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EDITORIAL 


The Third EuroConference on Nonlinear Dynamics in Physics and Related Sciences was 
held in Montecatini Terme, Tuscany, Italy, from 16-18 May, 1997, under the Auspices of a 
HCM granting programme (EEC Contract N. ERB4050PL934039). The title of the Confer¬ 
ence was: Control of Chaos: New Perspectives in Experimental and Theoretical Nonlinear 
Science. 

The aim of the workshop was to bring together researchers from different fields working 
on control and synchronization of chaos and spatiotemporal chaos, in order to assess the 
recent achievements in these areas, and to point out future directions. 

The subjects covered in the conference include the control and tracking of unstable 
periodic orbits in experimental systems, the synchronization of chaos for secure communi¬ 
cation systems, the control of chaotic laser dynamics, the stabilization of unstable periodic 
patterns in extended dynamical systems, the control of chemical reactions, and the targeting 
of chaotic trajectories to desired states. 

The control of chaos refers to a process wherein a chaotic dynamics is slightly perturbed 
in order to stabilize some of the unstable periodic orbits which are embedded within the 
chaotic attractor. This process allows us to use a single chaotic system to produce an 
infinite variety of periodic behaviors, with great flexibility in switching from one to another. 
Chaos synchronization refers to a situation where two identical chaotic systems starting 
from different initial conditions are kept synchronized by means of the transmission of a 
signal. Finally, targeting of chaos is a process wherein tiny perturbations applied to the 
system are able to steer the trajectory emerging from a given initial condition to any desired 
point of the attractor in a finite (and usually very short) time. 

In the last few years, these subjects have attracted more and more interest within the 
scientific community due to the extreme interdisciplinary nature of these research fields. For 
instance, applications of chaos synchronization range from building a secure communication 
system to the problem of bacterial replication. On the other hand, control of chaos has 
been experimentally tested in different dynamical situations, such as in laser dynamics, in 
electronic circuits, in the dynamics of the magnetoelastic ribbon, in the control of cardiac 
timing, in the control of chemical waves. 

All the different achievements have faced different aspects of the same problem, and in 
many cases different methods have been used to control different dynamical situations. 

More than 60 participants of the conference have profusely discussed these subjects 
from diverse points of view: from applied mathematics, to engineering, to laser physics, 
to chemistry, to electronics. We believe that the Montecatini Workshop has strongly con¬ 
tributed by bringing together different skills and expertises, and emphasizing the common 
achievements and the future directions of the field. 

The present special issue summarizes the discussion that took place during the confer¬ 
ence, by reporting most of the oral and poster presentations. This special issue also contains 
two tutorial papers, which are intended for the nonexpert readers as an introduction to the 
field. 

We gratefully thank all colleagues who helped referee the papers in this issue. Their 
comments, suggestions and remarks have significantly improved the quality of the presenta¬ 
tion of the reported papers. We would like to thank the EU representative, and all colleagues 
who chaired the different sections of the conference, who helped in putting up a successful 
and fruitful workshop. Special thanks also go to all the participants of the workshop and 
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to all the contributors to this special issue. We are particularly indebted to the following 
institutions for their contribution to the success of the conference: 

— United States Air Force European Office of Aerospace Research and Development 
(E.O.A.R.D.) 

— Istituto Nazionale di Ottica (Florence, Italy) 

— Gruppo Nazionale di Elettronica Quantistica e Plasmi del Consiglio Nazionale delle 
Ricerche (G.N.E.Q.P.) 


F. T. Arecchi 
S. Boccaletti 
M. Ciofini 
G. Grebogi 
R. Meucci 
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Tutorials and Reviews 


THE CONTROL OF CHAOS: THEORETICAL SCHEMES 
AND EXPERIMENTAL REALIZATIONS 

F. T. ARECCHF, S. BOCCALETTI*, M. CIOFINI and R. MEUCCI 
Istituto Nazionale di Ottica, 1-50125 Florence, Italy 
* Department of Physics, University of Firenze, 1-50125 Florence, Italy 
' Department of Physics and Applied Mathematics, 

Universidad de Navarra, E-31080, Pamplona, Spain 

C. GREBOGI 

University of Maryland, College Park, USA 
Received July 31, 1997; Revised November 28, 1997 


Controlling chaos is a process wherein an unstable periodic orbit embedded in a chaotic attractor 
is stabilized by means of tiny perturbations of the system. These perturbations imply goal 
oriented feedback techniques which act either on the state variables of the system or on the 
control parameters. We review some theoretical schemes and experimental implementations for 
the control of chaos. 


1. Introduction 

Controlling chaos consists in perturbing a chaotic 
system in order to stabilize a given unstable pe¬ 
riodic orbit (UPO) embedded in the chaotic at¬ 
tractor (CA). The UPO’s constitute the skeleton of 
chaotic dynamics, which, indeed, can be seen as a 
continuous irregular jumping process among neigh¬ 
borhoods of different periodic behaviors [Auerbach 
et al., 1987]. Thus, control of chaos implies the ex¬ 
traction of desired periodic motions out of a chaotic 
one, through the application of judiciously chosen 
small perturbations. The process allows to ex¬ 
ploit a single dynamical system for the production 
of a large number of different periodic behaviors, 
with an extreme flexibility in switching from one 
to another, so that the single system can carry out 
different performances with different yields. 

The aim of this paper is to summarize some the¬ 
oretical and experimental implementations of the 
above concepts. It is important to point out that 
the body of literature on this topic is very wide, 
and that the methodologies described here are by no 
means the only ones valid. Nevertheless, we hope 


that the reading of this paper may help researchers 
in entering this field, and in getting their bearings 
among the different methods. 

In Sec. 2 we report the first control method, 
which was proposed by Ott et al. [1990] (OGY) 
and which consists in slight readjustments of a 
control parameter each time the trajectory crosses 
the Poincare section (PS). Since a generic UPO is 
mapped on the PS by an ordered sequence of cross¬ 
ing points, OGY is able to stabilize such a sequence 
whenever the chaotic trajectory visits closely a 
neighborhood of one of the saddle PS points. The 
time lapse for a natural passage of the flow within 
the fixed neighborhood (hence for switching on the 
control process) can be very large. To minimize 
such a waiting time, a technique of targeting has 
been also introduced [Shinbrot et al, 1993]. 

Another technique to constrain a nonlinear sys¬ 
tem x(t) to follow a prescribed goal dynamics g (t) 
[Plapp & Huebler, 1990] is based upon the addition 
to the equation of motion dx/dt = F(x) of a term 
U(t) chosen in such a way that |x(f) — g(f)[ -4 0 
as t —> oo. Plapp and Huebler choose U(t) = 
dg/df — F(g(f)). The method provides robust 
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solutions, but in general the perturbation U would 
not be a small portion of the unperturbed dy¬ 
namics F. 

In other papers the effects of periodic [Lima 
&; Pettini, 1990; Braiman & Goldhirsch, 1991; 
Azevedo & Rezende, 1991] and stochastic [Fahy & 
Hamann, 1992] perturbations is explored in produc¬ 
ing dramatic changes in the dynamics, which even¬ 
tually may lead to the selection of some UPO, even 
though they cannot be considered in general as goal 
oriented. This parametric perturbation avenue is 
explored in several companion papers of this issue. 

Section 3 presents a method based upon the 
continuous application of a delayed feedback term 
in order to force the dynamical evolution of the sys¬ 
tem toward the desired periodic dynamics whenever 
the system visits such a periodic behavior closely 
[Pyragas, 1992]. 

Section 4 describes a frequency-domain control 
technique called “washout filter” [Tesi et al., 1996; 
Basso et al., 1997], based upon the insertion of a 
selective filter within a feedback loop. 

Section 5 introduces the method of adap¬ 
tive recognition [Arecchi et al., 1994] and control 
[Boccaletti & Arecchi, 1995, 1996] of chaos. The 
method has later been successfully applied to chaos 
synchronization [Boccaletti et al., 1997a], targeting 
of chaos [Boccaletti et al, 1997b], filtering noise 
from chaotic data sets [Boccaletti et al., 1997c], 
and eventually to the quenching of defects in an 
infinite dimensional dynamical system [Boccaletti 
et al, 1997d]. The technique consists of a first step 
wherein the unperturbed features of the dynamics 
are extracted, and in a second step in which pertur¬ 
bations are done for the control of desired periodic 
orbits. 

In Sec. 6 we review a few implementations of 
control of chaos in several experimental situations. 
The first experimental application of OGY was the 
stabilization of periodic orbits of a chaotic gravi¬ 
tationally buckled, amorphous magnetoelastic rib¬ 
bon [Ditto et al, 1990]. OGY inspired an eas¬ 
ily realizable experimental technique called OPF 
(occasional proportional feedback), and demon¬ 
strated in a chaotic diode oscillator [Hunt, 1991]. 

Many other experimental systems have pro¬ 
vided successful examples of chaos control. We re¬ 
call among the others the thermal convection loop 
[Singer et al, 1991], the yttrium iron garnet os¬ 
cillator [Azevedo & Rezende, 1991], the optical 
multimode solid-state laser [Roy et al, 1992], the 
Belouzov-Zabotinsky chemical reaction [Peng et al., 


1991; Petrov et al, 1993], the optical fiber laser 
[Bielawski et al, 1993], the CO 2 laser with modula¬ 
tion of losses [Meucci et al, 1994], and the lead-salt 
diode laser [Chin et al, 1996]. 

Finally, we discuss the experimental imple¬ 
mentations of the “washout filter” technique of 
Sec. 4. The method has provided successful experi¬ 
mental applications both on autonomous [Meucci 
et al., 1997] and nonautonomous [Ciofini et al, 
1997] chaotic CO 2 laser systems. 

2. The OGY Control of Chaos 

In this section we summarize the OGY method for 
the control of chaos [Ott et al, 1990]. Even though 
such method holds regardless of the number of pos¬ 
itive Liapunov exponents, for simplicity, we refer 
to a continuous-time chaotic dynamical system in 
a three-dimensional phase space, thus with a single 
positive Liapunov exponent. This is ruled by the 
differential equation 

f = F(x.p), (1) 

where x is a D dimensional vector (D ~ 3), and p is 
a control parameter that we assume to be accessible 
for adjustments. The goal is to temporally program 
such adjustments so as to achieve stabilization of 
some UPO embedded within the chaotic attractor. 
Furthermore, we imagine that the functional form 
of F is not known, but that experimental time se¬ 
ries of some scalar component z(t ) can be measured. 
By means of time-delay coordinates, and selecting 
an embedding time T, it is possible to reconstruct 
a M + 1 dimensional embedding space containing 
the vectors of the form 

X(<) = [z(t), z(t - T), z(t - 2 T), ..., z(t- MT )]. 

If one is interested in periodic orbits, one shall 
use X to obtain a surface of section, wherein 
any continuous-time-periodic orbit emerges as a 
discrete-time orbit cycling through a finite set of 
points. The requirement is that the embedding 
space has as many dimensions as there are coordi¬ 
nates of the point (M = D — 1), so that our surface 
of section is, in the present case, a two-dimensional 
surface. Let us now suppose that the control pa¬ 
rameter p can be varied in a small interval about 
some nominal value po (in the following, and with¬ 
out loss of generality, we take po = 0), ranging 
within the interval p* > p > — p*. Again, let us 
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suppose that the experimental measurement con¬ 
tains many points in the surface of section for p = 0. 
By denoting such points as £ 1 , £ 2 , £ 3 , • • •, £fc, then 
£ n are in general the coordinates at the nth inter¬ 
section of the surface by the orbit X(f). A common 
choice of the surface is z(t — MT) = constant, so 
that £„ = [z(t n ),..., z(t n - (M - 1)T)], t = t n de¬ 
noting the instant of the nth intersection. 

Finally, let us denote by £ = £/? = 0 a given 
fixed point, by X s and A u respectively the stable and 
unstable eigenvalues of the surface at £f, by e s and 
e u the unit vectors of the corresponding experimen¬ 
tally determined eigenvectors. When a change in p 
from p = 0 to some other value p = p is done, then 
the fixed point coordinate will shift correspondingly 
from 0 to some nearby point £f(|>)- For small p the 
following approximation holds 

<9£f (p) 

G ~ dp 

which allows an experimental estimate of G. On 
the surface and near £ = 0, we can describe the 
dynamics with the linear map 

£n+i - £f(p) - M • [£„ - £f(p)] 

where M is a 2 x 2 matrix. Using £f(p) — pG, the 
above equation reads 



£n+i — PnG [A„e u f„ + A s e s fs] • [£ n p n G] (2) 

where f s and f u are controvariant basis vectors de¬ 
fined by f s • e s = f„ • e„ = 1, f s • e u = f u • e s = 0. 
Let us assume £ n be located within a neighborhood 
of the desired fixed point. The control method con¬ 
sists in selecting p n so as £ n+ i be put onto the stable 
manifold of £ = 0, which implies to select p n so that 
f„ • £ n +i = 0. When £ n+ i falls on the stable mani¬ 
fold of the desired fixed point, the parameter can be 
set again to p = 0, because, by subsequent natural 
evolution, the dynamics will approach the desired 
fixed point at a geometrical rate X s . 

Dotting (2) with f u , we obtain 

_ ^ u fu ' £n 

Pn ~ X u - 1 T u • G 

which can be used provided the magnitude of the 
right-hand side be smaller than p*. In the oppo¬ 
site case, p n is set to 0. As a consequence, the 
perturbation is activated only when £ n is located 
within a narrow strip |£“| <£*,£“ = f u • £ n and 
£* = P*|(l - A~ 1 )G • f u |. 


This way, a stable periodic orbit is obtained 
out of the chaotic evolution of the dynamics. As 
we mentioned, the control of chaos gives flexibility. 
By turning the small controlling perturbations off, 
one can switch the time asymptotic behavior from 
one periodic orbit to another. In some situations, 
where flexibility offered by the ability to do such 
switching is desirable, it may be advantageous to 
design the system so that it is chaotic. In other sit¬ 
uations, where one is presented with a chaotic sys¬ 
tem, the method may allow one to eliminate chaos 
and achieve great improved behavior at relatively 
low cost. The OGY ideas can also be applied to 
stabilize a desired chaotic trajectory, which has po¬ 
tential applications to problems such as synchro¬ 
nization of chaotic systems [Lai & Grebogi, 1993], 
conversion of transient chaos into sustained chaos 
[Lai & Grebogi, 1994], communication with chaos 
[Hayes et al., 1993, 1994; Rosa et al ., 1997; Bollt 
et al., 1997], and selection of a desired phase [Nagai 
& Lai, 1995]. 


3. The Pyragas’ Method 

An alternative time-continuous method [Pyragas, 
1992] considers a dynamical system ruled by a set of 
unknown ordinary differential equations, and hav¬ 
ing some scalar variable accessible for measure¬ 
ments. Furthermore, the system possesses at least 
one input accessible for external forcing. The above 
assumptions are met by the following model 

^ = P(y, x) + F(t ); ^ = Q {y, x) (3) 

where y represents the output scalar variable, x the 
remaining hidden variables of the dynamical sys¬ 
tem, F(t) is an input signal which disturbs the dy¬ 
namical evolution of the variable y , and P and Q 
are two nonlinear functions. 

Let us imagine that system (3) produces chaotic 
dynamics for F = 0. In general, a large num¬ 
ber of the UPO’s within a chaotic attractor can 
be obtained from a single scalar variable through 
the standard method of delay coordinates. There¬ 
fore, one can extract from the measured variable 
y various periodic signals of the form y = yi{t), 
yi(t + Ti) = yi(t), where T t represents the period of 
the ith unstable periodic orbit. 

To achieve stabilization of the selected UPO, 
let us design an external feedback line which rein¬ 
jects into the system the difference D(t) between 
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the signals y(t ) and y(t — r) as a control signal: 

F(t) = K[y(t - t) - y(t)\ (4) 

where the weight K has to provide a negative feed¬ 
back (K < 0) and r represents a time delay. Stabi¬ 
lization of the ith UPO is achieved when r equals 
the period T). 

The method can be extended by adding infor¬ 
mation on previous periods, that is, replacing (4) 
with [Bleich k. Socolar, 1996] 

( OO 

fc=l 

( 5 ) 

with 0 < R < 1 a suitable real parameter, and k 
integer. 

Notice that use of the perturbations (4) or (5) 
transforms the Ordinary Differential Equation (3) 
to a Delay Differential Equation. This requires 
some warning. As a delayed dynamical system is 
richer than an instantaneous one, care should be 
put in stabilizing a true UPO of the original unper¬ 
turbed system, rather than a spurious UPO intro¬ 
duced artificially by the delay. 

4. The Washout Filter 

Control of chaos can be achieved by negative feed¬ 
back of suitable spectral components of a system 
variable [Meucci et al., 1996; Tesi et al., 1996; Basso 
et al., 1997; Ciofini et al., 1997]. The set up con¬ 
sists in a feedback loop wherein all unwanted fre¬ 
quencies present in the chaotic spectrum are trans¬ 
mitted as correction signals by means of a selective 
filter (“washout filter”). In this way, the system is 
allowed to oscillate at the only frequency which is 
not fed back, namely, that of the unstable orbit to 
be stabilized. This control scheme is very easy to be 
implemented, besides having the relevant advantage 
of being robust and fast. 

Consider a dynamical system modeled by ordi¬ 
nary differential equations in the form 

±i=F(xi,x 2 ) ^ 

±2 = G(x i, x 2 ), 

F and G being a linear and a nonlinear function, 
respectively. The forthcoming considerations still 
hold in case the system is nonautonomous. Consid¬ 
ering a stationary periodic regime where each vari¬ 
able can be approximated by Xk ~ e st (s = ius), let 
us introduce in the first equation a suitable nega¬ 
tive feedback loop for x\ through a “washout filter”, 



characterized by a certain transfer function C(s ): 
sxi = F(x i, x 2 ) — x\C(s). 


The above equation can be also rewritten as 
F(x i, x 2 ) 

* 1 = 7 * 


In order to stabilize a given orbit with pulsation 
Q introducing a minimal perturbation, C(s) should 
vanish for u = 0 and uj = Q (which implies that the 
feedback does not alter the fixed point and the limit 
cycle solutions of the unperturbed system). More¬ 
over, depending on the route to chaos, it is useful to 
choose C(s ) so that it presents a maximum in cor¬ 
respondence to the frequency characteristic of the 
transition to chaos. Whenever the above conditions 
are fulfilled, the only frequency component which is 
not affected by the feedback is that of the cycle to 
be stabilized, while all the other components are 
sent back as a correction signal. 

For example, in the case of the subharmonic 
route to chaos, the filter structure for stabilizing the 
period-1 (fundamental frequency / = orbit 

can be modeled as 


C(s) = (3- 


s(s 2 + D 2 ) 


s 2 + £Ds + j 


(s + 


where £ = 0.4, fi = 1.5 and (3 is the gain factor. The 
amplitude and phase responses of the above trans¬ 
fer function are shown in Fig. 1 for / = 110 kHz; 
the maximum of C(s ) is set approximately at U/2. 

We add some remarks on the applications to 
real experimental conditions. 


(i) This control scheme can be in principle very 
fast; indeed, the feedback loop can be entirely 
realized by analog electronics. 

(ii) The control is also very robust since it is inde¬ 
pendent of parameter fluctuations. 

(iii) Regarding the possibility of stabilizing more 
complex orbits (i.e. a period-2 orbit or a 
torus), one has just to design a different fil¬ 
ter C(s), with several zeroes corresponding to 
all the frequency components of the cycle to be 
stabilized. 

(iv) The basic structure of Eqs. (6) (called Lur’e 
form) is peculiar of dynamical systems widely 
studied in the literature, such as the Chua cir¬ 
cuit, the Rossler model and the Duffing oscil¬ 
lator [Genesio et al., 1993]. 
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Fig. 1. Amplitude and phase response curves of the washout filter as a function of oj. 


As a final consideration, we observe that 
this method presents several similarities with the 
time-delayed autosynchronization method [Pyra- 
gas, 1992; Basso et al., 1997], despite the conceptual 
difference between these two techniques. In the fre¬ 
quency domain, the Pyragas’ method can be seen 
as a negative feedback loop with a high-pass filter 
whose amplitude response goes to zero at the fun¬ 
damental frequency Cl and for all its multiples nCl. 
In this way, the stabilized limit cycle is exactly the 
same as that of the unperturbed system since it 
contains all the harmonics, while the perturbation 
vanishes. 


5. The Adaptive Algorithm 

Let us consider a dissipative dynamics ruled by 
Eq. (1). The adaptive control technique consists 
in two successive steps, the first one, in which the 
unperturbed features of the dynamics are extracted 
[Arecchi et a/., 1994] and the periods of the UPO’s 
are measured, and the second one whereby adaptive 
perturbations are applied in order to stabilize the 
selected UPO [Boccaletti &; Arecchi, 1995, 1996]. 

We consider an observer “blind” to the main co¬ 
ordinate position Xi (i = 1,..., D) and interested 
only in its variation 


Sxi(t n - |-i) — Xi(t n + 1) Xi(t n ) , (7) 

where t n+ i —t n = T n is the nth adjustable interval, 
to be specified. In order to assign T n +\ we consider 
the local variation rate 


Ai(tn+l) — In 
Tn 


$%i{tn+ 1 ) 

fixi (tn) 


( 8 ) 


hi 

Here r n+ i is the minimum of all corresponding 
to the different i, and defined by the rule 

T nli = T nH l ~ tanh(CTAj(t n +i))). (9) 

Equation (9) arises from the following consider¬ 
ations. To obtain a sequence of geometrically regu¬ 
lar Sxi , we shrink (stretch) the time intervals when¬ 
ever the actual value of 5xi is larger (smaller) than 
the previously observed one. The hyperbolic tan¬ 
gent function maps the whole range of <rX l into the 
interval (—1, +1). The constant a, strictly positive, 
is chosen in such a way as to forbid r -jp from going 
to zero. It may be taken as an a priori sensitivity, 
however, a more sensible assignment would consist 
in fixing cr by a moving average procedure, looking 
at the unbiased dynamical evolution for a while and 
then taking a cr value smaller than the reciprocal of 
the maximal A recorded in that time span. Notice 
that a moving sensitivity is mandatory whenever 
the adaptive recognition is specialized to the mea¬ 
surement of a periodic orbit [Arecchi et al ., 1994]. 

We thus obtain a sequence of observation times 
starting from to 

to,h = t 0 + f, t 2 = t\ + 7l, . .., t n+ j = t n + r n ,.. . 

( 10 ) 

corresponding to which the variations of 5xi(t n ) can 
be reduced below a preassigned value. 

The observations performed at these times pro¬ 
vide a “regularized” window, and the time sequence 
(10) now includes the chaotic information which was 
in the original geometric sequence x(i). The se¬ 
quence (10) contains the relevant information on 
the dynamics, and we can characterize chaos as 
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Fig. 2. Return maps r n +1 versus r n for (a) Ro4 and (b) Ro4 with an additional 1% white noise. Initial conditions: a:(0) = —20, 
y{ 0) = z(0) = 0, w(0) = 15. g = 0.000048. 


follows. Since in Eq. (9) |<rAj| <C 1, then two suc¬ 
cessive r n must be strongly correlated. As a result, 
even though the set of r n may be spread over a 
rather wide support, the return map r n+ i versus T n 
must cluster along the diagonal. The residual devi¬ 
ations, averaged over many successive n’s, provide 
the decorrelation trend, and hence yield an accu¬ 
rate assignment of the maximum Liapunov expo¬ 
nent. However, presence of large deviations from 
the diagonal denotes an uncorrelated perturbation. 
This may be some additive noise, which eventually 
can be filtered out, thus extracting the determinis¬ 
tic dynamics [Boccaletti et al. , 1997c]. 

In the following we will summarize the ap¬ 
plication of such a method to the Roessler four¬ 
dimensional (Ro4) model [Roessler, 1979] for a 
vector state x = (aq, aq, £ 3 , £ 4 )- For particular 
initial conditions (aq( 0 ) = — 20 , £ 2 ( 0 ) = 0 : 3 ( 0 ) = 0 , 
£ 4 ( 0 ) = 15) and control parameters, Ro4 under¬ 
goes a hyperchaotic dynamics with two positive 
Liapunov exponents. 

Figure 2 reports the return map of the r n for 
Ro4 and for Ro4 with 1% noise. 

Since we are interested in stabilizing periodic 
dynamics, we need to extract the periods of UPO’s 
embedded within the CA. For this purpose, instead 
of considering the single step map, we construct the 
maps r n+ fc versus r n , k = 1 , 2 ,... and we plot the 
r.m.s. rjik) of the point distribution around the di¬ 
agonal of such maps as functions of the step interval 
k. For chaotic dynamics, temporal selfcorrelation 
lasts only for a finite time, hence one should ex¬ 
pect to obtain a monotonically increasing function 


rj(k). In fact, the chaotic dynamics steers the phase- 
space trajectory toward neighborhoods of different 
UPO’s. As the trajectory gets close to an UPO of 
period Tj, temporal selfcorrelation is rebuilt after 
Tj and the distribution of r includes windows of 
correlated values appearing as minima of 77 versus 
k around kj = Tj / (t) , (r) being the average of the 
r distribution. 

To give an exemple, we report in Fig. 3 the 
77-fc plot for Ro4, from which one can extract the 
different UPO’s periods looking to the minima of 
the 77 curve. In fact, during the observation, the 



Fig. 3. 77 -fc plot for Ro4 attractor. Initial conditions and 

control parameters as in the text. The recognition task has 
been performed with g = 0.01. Vertical axis has to be multi¬ 
plied for 10 -4 . 
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time intervals are changing, and a more rigorous 
determination of the period is provided by mini¬ 
mizing a suitable cost function in the vicinity of 
each minimum of the 77 curve [Boccaletti & Arecchi, 
1995]. 

Once the periods Tj (j = 1 , 2,...) of the UPO’s 
have been measured, stabilization of each one can 
be achieved when the system naturally visits closely 
phase space neighborhoods of that UPO. For a 
nonauthonomous system, a period T may corre¬ 
spond to many degenerate UPO’s. In this case, se¬ 
lection of the desired one can be achieved by the 
study of the topology of all UPO’s corresponding 
to the same period and by switching-on the con¬ 
trol task when the system is shadowing the selected 
UPO. Several topological approaches to the UPO’s 
detection have been developed [Cvitanovic et al ., 
1988; Gunaratne et al., 1989; Mindlin et al., 1990; 
Tufillaro et al., 1990]. 

The control procedure is done by use of the 
following modified algorithm. At each new ob¬ 
servation time t n+ 1 = t n + T n and for each com¬ 
ponent i of the dynamics, instead of Eq. (7), we 
evaluate the differences between actual and desired 
values 

5Xi(t n +l) = Xj(t n +x) ~ x i{tn- 1-1 — Tj) , (11) 

and the local variation rates A’s now read 


Ai(£ n+ i) = — log 

T n 


x i{tn+ 1) x i{tn +1 Tj) 
x i{ln) x i{tn Tj) 


( 12 ) 


We keep Eq. (9) and the choice of the minimum 
for the updating process of r’s. Defining U(f) as 
the vector with ith component (constant over each 
observation time interval) given by 

Ui(t n + 1 ) = (xj(t n _i_i — Tj) — Xj(f n -f 1 )), (13) 

x n+l 

we add such a vector to the evolution equation, 
which then becomes 

5|=F(z,p) + U((). (14) 

Notice that the A’s given by Eq. (12) track the 
rate of separation of the actual trajectory from the 
desired one. Indeed, A negative means that locally 
the true orbit is collapsing into the desired one and 
hence the actual dynamics is shadowing the desired 
UPO, while A positive implies that the actual tra¬ 
jectory is locally diverging away from the desired 
one and control has to be performed in order to 
constrain the orbit to shadow the desired UPO. 

As a consequence, contraction or expansion of 
r’s now reflects the need to perturb the dynamics 
more or less robustly in order to stabilize the desired 
UPO. This appears as a weight to the correction of 
Eq. (13), which, once a given Tj has been chosen 
by the operator, is selected by the same adaptive 
dynamics. 

Once again, the introduced adaptive weighting 
procedure in Eq. (13) assures the effectiveness of the 



Fig. 4. (a) (an, £ 2 ) projection of the phase space portrait for the controlled period-8 of Ro4 attractor. Control task has been 

performed with period-8 extracted from Fig. 3 and g = 10 -5 . (b) Time evolution of the first component x\ of Ro4 before and 
after control. Arrows indicate the instant at which control task begins. 
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time 

Fig. 5. (a) Temporal evolution of the first component of the additive controlling term U\ during the control of period-8 of 

Ro4 and (b) temporal evolution of the uncontrolled dxi/dt. The adaptive correction term is between two and three orders 
of magnitude smaller than the natural evolution of the dynamics. Same stipulation for controlling task as mentioned in the 
caption of Fig. 4. 


method (perturbation is larger or smaller whenever 
it has to be) as well as the fact that the additive 
term U is much smaller than the unperturbed dy¬ 
namics F. 

In Fig. 4 we show the control of period-8 of Ro4. 

Figure 5 reports the perturbation U\(t) and the 
unperturbed dynamics Gi for the Ro4 model dur¬ 
ing the control task of period-8, in order to show 
that the former is between two and three orders of 
magnitude smaller than the latter as expected by 
the above discussion. 

Notice that the limit a = 0 of the above al¬ 
gorithm recovers [Pyragas, 1992]. Choosing cr ^ 
0 implies an adaptive nature of the forcing term 
[Eq. (13)] which is inversely proportional to the time 
intervals and hence is weighted by the information 
extracted from the dynamics itself. 

Applications of this method have been already 
reported in the introductory section. 


6. Review of Some Control 

Experiments 

The OGY method described in Sec. 2 has found ap¬ 
plications in several experimental situations. As an 
example, we herewith report the first experimental 
realization of chaos control which was done in 1990 
by Ditto et al. [1990]. 

In this case, the theoretical background of 
Sec. 2 has been applied to an experimental system 
consisting of a gravitationally buckled, amorphous 
magnetoelastic ribbon. The ribbon exhibits large 


reversible . changes of its Young’s modulus when 
small magnetic fields are applied. It is clamped at 
the base, yielding a free vertical length greater than 
the Euler buckling length, thus providing an initial 
buckling configuration. The ribbon is placed within 
three mutually orthogonal pairs of Helmoltz coils, 
allowing to compensate for the Earth’s magnetic 
field. Finally, a uniform vertical magnetic field is 
applied along the ribbon. In this configuration, the 
Young’s modulus of the ribbon is varied due to the 
applied vertical magnetic field which has the form 
H — H dc + H ac cos(2nft). Both amplitudes have 
been set typically in the range 0.1-2.5 Oe. A mea¬ 
surement of the curvature of the ribbon near the 
base is provided by a sensor. 

The experimental data consist in time series 
voltages V (t) acquired from the output of the sensor 
and sampled at the drive period of the ac magnetic 
field (i.e. at times t n = n //). The sampled voltages 
are considered as iterates of the map X n = V(t n ) 
and the control theory in Sec. 2 is applied, taking as 
control parameter the amplitude of the continuous 
component of the magnetic field H dc . Selecting H ac , 
H dc and / so as to produce chaotic dynamics, con¬ 
trol of period-1 UPO is achieved for over 200 000 
iterates (approximately 64 hours) with maximum 
perturbation of about 1% of the unperturbed con¬ 
trol parameter. With the same setup, control of one 
of the period-2 UPO is also achieved. 

A modification of the OGY method, called oc¬ 
casional proportional feedback (OPF) has been used 
to stabilize unstable orbits in a chemical system 
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[Peng et al., 1991, Petrov et al., 1993] and in a diode 
resonator [Hunt, 1991]. This technique consists 
in feeding back deviations of the chaotic variable 
within a specified window from a reference point to 
perturb a control parameter. The chemical experi¬ 
ment deals with the Belousov-Zhabotinsky reaction 
carried out in a continuous-flow stirred-tank reac¬ 
tor. The flow rate p of the reactants into the tank 
ultimately determines whether the system shows 
steady state, periodic and chaotic behaviors. The 
control algorithm takes advantage of the predictable 
evolution of the chaotic system in the vicinity of a 
fixed point in the next-amplitude return map of a 
suitable variable A (dimensionless concentration). 
The position of the period-1 fixed point is given by 
the intersection of the map with the bisectrix, where 
A n+ 1 = A n = A s . The map can be shifted to tar¬ 
get the fixed point by applying a perturbation to the 
bifurcation parameter pi according to the difference 
between the system state and the fixed point 


where g is a suitable constant. In an analogous 
manner (changing A s and g) period-2 and period-4 
unstable orbits have been stabilized. 

The electronic experiment consists of a p-n 
junction rectifier in series with an inductor. The 
system exhibits the period-doubling route to chaos 
when driven with an increasing sinusoidal voltage. 
The current through the diode provides a conve¬ 
nient chaotic variable; if the peak current I n falls 
within a given window, the driving voltage is am¬ 
plitude modulated with a signal proportional to the 
difference between I n and the center of the window. 
By changing the level and the width of the window, 
or the gain of the feedback signal, several unstable 
orbits are stabilized, up to the period-23. 

The OPF has been also applied by Roy et al. 
[1992] to an autonomously chaotic multimode laser, 
that is, a high dimensional system for which the 
chaotic attractor is not characterized by a low¬ 
dimensional map. The experimental setup con¬ 
sists of a diode-laser-pumped solid state Nd-doped 
yttrium aluminum garnet (Nd:YAG) laser contain¬ 
ing a KTP doubling crystal. The source of chaotic 
behavior in this laser is the coupling of the longitu¬ 
dinal modes through the nonlinear process of sum- 
frequency generation. This process destabilizes the 
relaxation oscillations which are normally damped 
in a system without the intracavity KTP crystal. 
The total laser output intensity is sampled within 


a window of selected offset and width, and a signal 
proportional to the deviation from the center of the 
window is applied to perturb the injection current of 
the pumping laser diodes. With this configuration 
period-1, -4 and -9 limit cycles have been success¬ 
fully stabilized with relative perturbation amplitude 
less than 10%. 

The control scheme proposed by Pyragas 
(discussed in Sec. 3) has been first experimentally 
demonstrated in an electronic circuit [Gauthier 
et al., 1994] and in a modulated laser [Bielawski 
et al., 1994], In the first experiment, the setup is 
similar to that reported by Hunt, the only differ¬ 
ence being that the resonator has been modified to 
operate at higher frequency (10 MHz). The con¬ 
trol is derived by directing half of the voltage sig¬ 
nal (proportional to the resonator current) directly 
into one input of a summing amplifier, while the 
other half, delayed and inverted, is sent to the sec¬ 
ond input. The delay line consists of a cable with 
length precisely adjusted so that it provides a delay 
t corresponding to the period of the desired UPO. 
The control signal is reinjected into the resonator as 
an additive perturbation of the driving voltage. In 
this way, a close reproduction of the control Eq. (4) 
is achieved, allowing stabilization and tracking of 
period-1, -2 and -4 unstable limit cycles. 

The second experiment deals with a CO 2 laser 
with cavity loss modulation, obtained by driving 
an intracavity electro-optic crystal with an external 
sinusoidal voltage. In this case, the chaotic variable 
is the infrared (10 pm) laser light, monitored by a 
fast photovoltaic detector. The detector voltage is 
used to modulate a laser diode emitting at 845 nm 
so that a time delayed voltage can be obtained by 
propagating the laser diode light in a long fiber and 
detecting it at the end. Finally, the difference be¬ 
tween the CO 2 laser intensity and its delayed ver¬ 
sion is added to the modulation signal after suitable 
amplification. With such a configuration the unsta¬ 
ble period-1 orbit is stabilized and tracked along a 
wide range in the bifurcation diagram. 

The experimental implementation of the con¬ 
trol scheme based on the washout filter has been 
tested in the chaotic regimes of both a nonau- 
tonomous system [Meucci et al, 1996; Ciofini et al., 
1997] (a CO 2 laser with externally modulated 
losses) and of an autonomous system [Meucci et al., 
1997] (a CO 2 laser with intensity feedback), obtain¬ 
ing stabilization and tracking of different unstable 
periodic orbits, with perturbation amplitudes of the 
order of few percent. 
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The experimental setup consists of a single 
mode CO 2 laser with an intracavity electro-optic 
crystal which can be driven by an external voltage 
V{t) to modulate the cavity losses k. In the nonau- 
tonomous case, the driving voltage V(t ) is a sinu¬ 
soidal signal, so that the cavity loss parameter k be¬ 
comes: k = ko[l + m sin(27r/f)], where / = 110 kHz 
is the modulation frequency and the modulation 
depth m, proportional to the amplitude of V(t), 
is the control parameter. By increasing m the sys¬ 
tem undergoes the transition to chaos through a se¬ 
quence of subharmonic bifurcations. The period-1 
orbit is stable up to m = 0.1, and a further increase 
of m drives the system to period-2 and period- 
4 orbits, followed by the first chaotic region and, 
finally, by a period-3 stable solution. The con¬ 
trol was implemented with a negative feedback loop 
where the laser intensity, revealed by a fast detector, 
is first filtered and then subtracted from V ( t ). The 


characteristic curves in Fig. 1 have been closely re¬ 
produced by a two-stage passive filter (a band-pass 
stage combined with a Notch stage) entirely realized 
by analog circuitry (Fig. 6). 

Figure 7 reports the experiment. Figure 7(a) is 
the chaotic laser oscillations (m = 0.18) observed 


*-2 



Fig. 6. Electronic scheme of the washout filter: Ri = 1 kfi, 
L\ = 6.5 mH, Ci = 0.1 nF, R 2 = 6.7 kfi, L 2 = 12.4 mH and 
C 2 = 0.2 nF. 



Time (jjls ) Time (ji s) 



Time (jjl s) Time (jjls) 


Fig. 7. Experimental results: (a) chaotic laser intensity without control and (b) corresponding control signal; (c) and (d) 
represent the same signals as (a) and (b), respectively, but in the case of activated control. 
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with open control loop, while Fig. 7(b) is the cor¬ 
responding control signal. The same signals are re¬ 
ported in Figs. 7(c) and 7(d), respectively, but in 
the case of a closed control loop. The control signal 
amplitude of Fig. 7(b) provides a 1.25% perturba¬ 
tion of the driving signal. This result confirms that 
the method allows stabilization of unstable orbits 
slightly different from those embedded in the un¬ 
perturbed system. 

Following the above criteria, a second filter has 
been prepared allowing stabilization of the period-2 
limit cycle, and both filters have been used to track 
the corresponding trajectory along the whole bifur¬ 
cation diagram. The unperturbed bifurcation dia¬ 
gram has been measured by increasing the control 
parameter m at fixed steps, and processing the laser 
intensity in order to extract the maxima (Fig. 8). 
The same measurements have been repeated after 
the insertion of the feedback loops with the two fil¬ 
ters. Figures 8(a) and 8(b) show the superposition 
of the unperturbed bifurcation diagram (dots) with 
the tracked period-1 and period-2 orbits (circles), 
respectively. In both cases the tracking has been 
achieved without any readjustment of the gain of 
the feedback loop over the whole explored range, 
and with relative perturbation amplitudes less 
than 3%. 

In a different experiment, the control has been 
tested on a CO2 laser made chaotic by an intensity 
feedback. In this case the system is autonomous, 
since the modulator voltage V(t) carries informa¬ 
tion on the output intensity. Indeed it is obtained 
by detecting the laser output and then amplifying 
such a signal. The equation for V ( t ) is 

where (3 is the damping rate of the feedback loop, 
I is the laser intensity, B a bias voltage acting as 
the control parameter and R the total gain of the 
loop (the term al represents the saturation of the 
detection apparatus). 

Starting from constant laser output and in¬ 
creasing B , the system passes to a limit cycle 
through a Hopf bifurcation, and then it reaches the 
chaotic region after a subharmonic cascade. The 
spectral analysis of the chaotic signal for B — 360 V 
shows the presence of a peak at / = 22 kHz, rem¬ 
nant of the Hopf bifurcation. This property sugr 
gests to prepare a suitable washout filter with zero 
amplitude in correspondence of /. Figure 9 shows 
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Fig. 8. Experimental results of the tracking of (a) period-1 
unstable orbit and (b) period-2 unstable orbit. Circles and 
points: maxima in the laser output signal with and without 
the control loop, respectively. 


the three-dimensional reconstruction of the chaotic 
attractor and the stabilized period-1 orbit obtained 
with a relative perturbation of about 7%. 

All the experimental results can be adequately 
reproduced by numerical integration of a CO2 laser 
model based on rate equations for the intensity and 
for the populations of the resonant levels coupled 
by collisions with the rotational manifolds. 
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Fig. 9. (a) Three-dimensional reconstruction of the chaotic 

attractor and (b) stabilized period-1 orbit. 
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In this paper the technical details of chaotic digital code-division multiple access ((CD) 2 MA) 
communication systems used in cable communication systems are presented. The cable com¬ 
munication system may be a pure coaxial RF cable network, a hybrid fiber-coax network, or 
a pure optical fiber network for high-capacity data link. As an example of its many poten¬ 
tial applications in cable communication systems, (CD) 2 MA is used to support the upstream 
digital data communications in cable TV systems occupying the very noisy 5-40 MHz portion 
of spectrum. Although the (CD) 2 MA proposed in this paper is only used to support current 
Internet services via cable TV networks, it can also be used to support the high-speed data-link 
provided in the 550-750 MHz band in hybrid fiber-coax networks. 

(CD) 2 MA is a new communication framework which uses band-limited chaotic carriers in¬ 
stead of linear ones. For the purpose of generality, in this paper the band-limited chaotic carriers 
are approximated by groups of linear sub-carriers, which distribute within the same bandwidth 
with a fixed amplitude, random phases and uniformly distributed frequencies. The theoretical 
result of the performance of (CD) 2 MA is given. We also provide the simulation results of the 
bit-error-rate (BER) performances of a synchronous (CD) 2 MA used in cable communication 
systems. The results show that the (CD) 2 MA system has a better performance than the syn¬ 
chronous CDMA system proposed for the same cable communication system. Technical details 
of (CD) 2 MA are also presented for the future design of prototype systems. We present the 
framework of the whole (CD) 2 MA system including carrier synchronization, timing recovery 
and the details of nonlinear carrier generators. 


1. Introduction 

Currently, the main networks connecting American 
homes to the Internet are telephone lines. However, 
the twisted pair copper telephone wire cannot pro¬ 
vide much more bandwidth than 56 to 64 kbps. For 
this reason, many American families have to install 
separate telephone lines for their PCs. 

Cable television system is a kind of communica¬ 
tion system which was originally designed to broad¬ 
cast television signals via coaxial RF cables rather 
than through the air [Baer, 1974], More than 50 


years older, the cable television system has covered 
almost every corner of North America and Europe. 
Today, cable TV companies have direct access to 
more than 63 million U.S. homes. Recently, due to 
the rapid growth of the Internet market, the func¬ 
tions of cable television systems have been changed 
from sending only analog television signals to send¬ 
ing both analog television signals and digital Inter¬ 
net information. In view of this conceptual change 
of usage of the cable television system, a brand new 
television set called WebTV (shown in Fig. 1) will 
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easy access with a 
remote control 


Fig. 1. The outline of a WebTV system. 


usher in a new era of television [Tomari et al, 1997], 
where entertainment and information services are 
combined into single devices. The main advantage 
of a cable TV network is its ability to transmit high- 
bandwidth video, voice and data. Furthermore, a 
cable TV network integrated with digital informa¬ 
tion sources can also create new Internet services 
such as interactive TV programs, high-speed on-line 
services and videophone services. 

Being a mature technical and widespread com¬ 
mercial service, it has become very expensive to up¬ 


grade the existing cable TV system. For example, 
by using optical fibers, the bandwidth of a cable TV 
system can be enlarged by 50%, from 40-550 MHz 
to 40-750 MHz. In addition, the optical fiber is 
more reliable than a coaxial' cable. However, even 
for a small city, it would cost $20 million to upgrade 
its cable system to fiber-optic lines. 

A less costly approach is to keep the hardware 
framework of the cable TV systems unchanged but 
exploit optimally the channel capacity of existing 
networks. One of these methods is synchronous 
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code division multiple access (S-CDMA) proposed 
by Terayon Communication Systems. 1 Terayon’s 
S-CDMA can provide a very robust transmission 
with a full 10 Mbps per 6 MHz channel upstream 
and downstream over a fiber-coax cable system. So 
far, since only a very small portion of the band¬ 
width, namely 5-40 MHz, can be used to provide 
two-way, high-speed data links between a subscriber 
and the Internet, the channel capacity of this por¬ 
tion of spectrum becomes very critical to the service 
quality for cable subscribers to access the Internet 
via cable TV networks. 

Besides upgrading the cable TV network, a ca¬ 
ble modem for converting the data stream to radio 
frequency should be plugged into a PC. The cur¬ 
rent price of a cable modem is about $250 to $600 
[Huffaker, 1995]. Since today’s PCs usually include 
pre-installed modems for accessing the Internet via 
telephone lines, a user may object to the extra cost 
of installing the cable modem. To overcome this ob¬ 
jection, a cable modem should be cheap enough for 
the cable TV company to provide each user with 
a.free cable modem at the outset. In designing 
the application of (CD) 2 MA to cable TV Internet 
services, we must always bear this consideration in 
mind such that all costly and sophisticated devices 
are concentrated at the fiber node, or at the head- 
end, to reduce the expense of the overall network. 

We have shown in [Yang & Chua, 1997] that an 
asynchronous (CD) 2 MA technology can double the 
channel capacity of an asynchronous CDMA system 
in a wireless communication environment. In this 
paper, we propose a synchronous (CD) 2 MA tech¬ 
nology and show that better performance than the 
S-CDMA can be achieved. Instead of employing 
a linear carrier, a (CD) 2 MA system uses a chaotic 
carrier. Whenever a chaotic carrier is used the non¬ 
linearity of the channel can be exploited to make 
the carriers more distinguishable from each other 
by reducing the correlation of the different chaotic 
carriers generated by the same chaotic generator 
structure in different transmitters. For linear com¬ 
munication systems, however, any nonlinearity will 
change the waveforms of the linear carriers such 
that it is more difficult to recover the modulated 
information because the current design principle 
is to maintain the waveform unchanged all along 
the channel. Although beautiful theoretical analy¬ 
sis can be formulated in view of the simplicity of 
linearity, the price a linear communication system 
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must pay is to ensure that all devices along the 
channel are operating only in linear regions. This 
makes the whole communication system expensive 
because nonlinearity is wasted as useless excess 
baggage. 

The organization of this paper is as follows. 
In Sec. 2, the layout of today’s cable TV system 
with Internet service is presented. In Sec. 3, The 
layout of the (CD) 2 MA system used in upstream 
and downstream cable TV systems is presented. 
In Sec. 4, the theoretical model of our proposed 
(CD) 2 MA scheme for cable TV systems is given. In 
Sec. 5, theoretical results on the performance of our 
(CD) 2 MA system are presented. In Sec. 6, a com¬ 
parison between the performance of the S-CDMA 
system and the (CD) 2 MA system for the 5-40 MHz 
band in the same cable TV system is presented. Fi¬ 
nally, some concluding remarks are given in Sec. 7. 

2. Cable TV Systems 

Figure 2(a) shows the outline of a typical current ca¬ 
ble communication system, which is not only a two- 
way analogue TV broadcast system but a two-way 
information network. However, unlike telephone 
networks, a cable TV network is very asymmetric, 
with over 90% capacity used for downstream and 
the rest for the upstream. The functional change of 
cable communication systems from one-way to two- 
way is driven by the rapid growth of the demand for 
bandwidth on the Internet. The central node of a 
cable TV system is the headend. Signals from differ¬ 
ent sources, such as satellite and terrestrial broad¬ 
cast, Internet, as well as local originating program¬ 
ming, are modulated onto radio frequency carriers 
and combined together for distribution over the ca¬ 
ble system. Supertrunks (high-quality microwave, 
fiber optic, or cable links) connect the head-end to 
local distribution centers, known as hubs. Several 
trunks may originate from a hub to provide cov¬ 
erage over a large contiguous area. Figure 2(b) 
shows a single trunk of a typical cable TV system. 
Tfunk amplifiers are installed along the trunk to 
maintain the signal level and compensate for ca¬ 
ble transmission characteristics. The bridge ampli¬ 
fier serves as a high-quality tap, providing connec¬ 
tion between the main trunk and multiple high-level 
branches. The line extender is a type of amplifier 
that maintains the signal level along the branch. 



splitter H 


line i 
extender 


tap tap sub-branch 


subscriber subscriber 


Fig. 2. The outline of typical current cable communication systems, (a) Block diagram of a typical cable communication 
system, (b) The block diagram of a single trunk. 
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The splitter connects a subbranch to a branch. Sub¬ 
scriber drops are connected to passive taps along 
the subbranch. Many cable TV companies support 
500 subscribers per fiber node in hybrid fiber-coax 
networks [HufFaker, 1995]. If each subscriber has 
a 64 kbps modem for its PC, then a 32 Mbps per 
fiber node should be designed for peak hours. How¬ 
ever, the available bandwidth is always swamped by 
more advanced Internet services such as teleconfer¬ 
ence where real time audio and video data streams 
are sent. 

Technical limitations on channel capacity are 
set principally by cable amplifiers. The coaxial ca¬ 
ble service occupies the 40-550 MHz portion of the 
spectrum, while the hybrid fiber-coax cable service 
occupies the 40-750 MHz portion of the spectrum. 
Since each commercial analog TV channel occu¬ 
pies a 6 MHz bandwidth, a cable TV system may 
provide up to 100 TV channels. In the future, a 
fully upgraded hybrid optical-coax cable network 
may use the 550-750 MHz band to provide digi¬ 
tal video, high-speed data and telephone services. 
In this spectrum band, (CD) 2 MA can also find 


applications. For a two-way cable TV system, the 
spectrum slot located at 5-40 MHz provides sub¬ 
scribers with an upstream data-link to the outside 
world through the fiber-coax cable system. For¬ 
mally, this upstream link can only provide a low 
speed data-link due to two main noise sources in 
this spectral slot (shown in Fig. 3). The first is 
impulsive noise from the electrical devices in PCs, 
TV sets, hair dryers, vacuum cleaner and etc., as 
well as vehicle emission systems. The second is 
narrowband interference picked up by the cable net¬ 
work itself such as Ham radio and Voice of America 
broadcasts. 

In the Terayon S-CDMA system, it was re¬ 
ported that a total capacity of 10 Mbps per 6 MHz 
channel can be achieved in the 5-40 MHz band. 
While this represents a big advance in technology, 
this channel capacity may be too small to match 
the rapid future growth of demand of the Inter¬ 
net service, accompanied by the occupation of PCs 
and WebTVs in subscriber homes. Although the 
soft-degradation of the service quality of current 
CDMA systems can provide a barely satisfactory 
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service to each subscriber during peak periods, such 
degradation of service can still be felt by subscribers 
if they are in the process of downloading image 
files, or some other data-intensive services, such as 
teleconferencing. 

3. (CD) 2 MA for 

Cable Communication System 

The main problems in using the 5-40 MHz band 
in cable TV systems are impulsive and narrowband 
noises. The multipath fading problem in wireless 
environments [Viterbi, 1995] does not occur in cable 
systems and the delay from each link can be mea¬ 
sured a priori Since the price that a subscriber is 
willing to pay for updating his cable service is rel¬ 
atively low, the transmitter and the receiver of the 
(CD) 2 MA system at the subscriber’s end should be 
as simple as possible. The layout of the (CD) 2 MA 
transmitter for upstream communication is shown 
in Fig. 4. Since the hardware structure of a ca¬ 
ble communication system is fixed during its opera¬ 
tion, each head-end controller can broadcast a local 
clock signal to synchronize the local clocks of all 
subscribers under the same headend. The timing 
recovery in a cable communication system is thus 
solved. 

The chaotic carrier used in the (CD) 2 MA sys¬ 
tem can be generated by an array of chaotic digital 
circuits, such as a “chaotic” (pseudo-chaotic) cel¬ 
lular automata [Toffoli, 1987], or a reversible cel¬ 
lular neural network [Crounse et al., 1996; Yang 
et al., 1996]. For manufacturing convenience, the 


hardware for each chaotic digital circuit is chosen 
to be the same in every transmitter. To ensure that 
at every moment the chaotic carriers in the same 
channel are orthogonal to each other, the headend 
controller must dynamically assign each transmitter 
a user ID number as the initial states of the digi¬ 
tal chaotic generator whenever a subscriber begins 
sending information. Then, the output of the dig¬ 
ital chaotic generator (which may be 8 or 12 bits) 
is fed into a D/A converter, whose output is sent 
to a frequency step-up for transferring the chaotic 
carrier to a different spectral band. Since the D/A 
converter may work at a clock speed of 1 MHz, we 
may have to shift this output to some prescribed 
6 MHz channel within the 5-40 MHz band. 

A frequency step-up is used to transfer the non¬ 
linear carrier to an appropriate spectral band. This 
is accomplished by a frequency multiplier , which 
consists of a nonlinear circuit followed by a band¬ 
pass filter, as shown in Fig. 5(a). While there are 
many choices for the nonlinear device, the simplest 
one is a transistor biased in the nonlinear region. 
The frequency multiplier based on a transistor is 
shown in Fig. 5(b). If a bandpass signal v\ n (t) is fed 
into a frequency multiplier, the output v out (t) will 
appear in a frequency band at the nth harmonic 
of the input frequency(range). However, the mul¬ 
tiplication factor n that can be provided by this 
circuit is usually small with a typical value of 3 
or 5 because the nonlinearity of a transistor is too 
“smooth”. To get a large gain, we need to find some 
device which has much more irregular nonlinearities 
such as breakpoints. 


user 

ID# 





Fig. 4. The block diagram of the upstream transmitter at subscriber’s end. 
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nonlinear 

device 



frequency multiplier 


v out(t) 


(a) 



(b) 

Fig. 5. A frequency multiplier used for the frequency step- 
up. (a) The block diagram of a frequency multiplier, (b) The 
circuit implementation of a frequency multiplier. 


The block diagram of the receiver at the head- 
end controller for retrieving the digital signal from 
the ith user is shown in Fig. 6. Since the delay of 
the ith subscriber has been previously measured, 
the receiver in the headend controller can find the 
pre-measured time delay from a lookup table and 
generate the corresponding chaotic carrier with the 
corresponding delay. The signal received by the ca¬ 
ble network is a mixture of noises and interferences 
from the other subscribers. The received signal is 
then multiplied by the regenerated carrier at the 
receiver. The result is low-pass filtered by a [0, T] 
integrator and then thresholded and sampled to give 
the recovered digital signal. The timing signal used 
in sampling is provided by the local clock signal 
at the head-end controller. This recovered signal 
is then sent out to the Internet from the headend 
controller. 

The downstream communication is almost the 
same as that of the upstream except that before the 
headend sends a chaotic carrier to user i, a time de¬ 
lay is compensated at the transmitter end. By doing 
this, the receiver at the subscriber’s end is almost 
the same as its transmitter except for an additional 
[0, T] integrator, a thresholding and a sampling 
circuit. The block digram of the receiver at the 
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Fig. 6. The block diagram of the receiver at the headend controller for upstream communication. 
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Fig. 7. The block diagram of the receiver at the subscriber’s end for downstream communication. 


subscriber’s end is shown in Fig. 7. Observe that 
the main part of the transmitter is shared by the 
receiver. Since the synchronization signal is 
broadcast through the whole network, the expen¬ 
sive synchronization circuit is not installed in the 
subscriber’s transmitter/receiver. Hence, the addi¬ 
tional equipment in the subscriber’s home is rela¬ 
tively inexpensive. 

Under ideal conditions, the power levels of the 
chaotic carriers from different users at the fiber 
node (or the headend) should be the same. This 
is the power control problem, which had been ex¬ 
tensively studied in wireless CDMA systems [Lee, 
1989]. However, in view of the fixed hardware 
structures of cable TV networks, amplifiers can be 
designed to satisfy this condition. In this paper, 
perfect power control is assumed. 

4. Model of Chaotic Digital CDMA 

Communication Systems 

The difference between classical CDMA and 
(CD) 2 MA is that the former uses a pseudo-random 
digital NRZ signal as a nonlinear sub-carrier and 
this carrier is then used to modulate a linear carrier 
for central frequency shift. The latter uses chaotic 
carriers directly in the desired spectral band. To 
clarify the difference between them, in this section 
we present examples to show the different principles 
used in these two schemes. 

Since in (CD) 2 MA systems many candidates of 
chaotic RF waveforms can be used as chaotic carri¬ 


ers, we give a theoretical model for a chaotic carrier 
which has a flat enough spectrum in the bandwidth 
we are interested in. The chaotic carrier for the ith 
user is given by 

N 

Ci(t) = Yl H n(f) COS(0 ) n t + Oi, n) 

71=1 

= Re jxj 0i > n (*) ex p(M*<+j0i,n) j , N^oo 

(1) 

where j = \/— 1, and a^ n (t) G 3ft, ui n G [prescribed 
spectral band] and G ( — 7r, 7r] are respectively 
the amplitude, frequency and phase of the nth com¬ 
ponent of the chaotic carrier. Since Ci(t) is a non¬ 
linear wave, at least one of a» >n , n = 1,..., N, must 
be time-varying. Observe that the single tune car¬ 
rier is a special case of c*(t) with an ^ 0 and = 0 
ifj^l for a given l. 

Similarly, the carrier for the jth user is modeled 

by 

N 

c j{t) = ^2 a j>n (t) cos(w n f + 9j t n) 

71=1 

N 1 

X) a j,n(t) exp(jUnt+jdj !n ) > , iV->00 

n= 1 J 

( 2 ) 

Since in (CD) 2 MA systems, we usually use a [0, T] 
integrator as the low-pass filter (LPF), the moving 
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average of the cross-correlation between a(t) and Cj(t), i / j with time window T is given by 
1 f l 

nj(t ) = - / Ci(t)cj(t)dt 

i Jt-T 

l ft N N 

= Pft I ^ ' ®j >n (f) COs(u) n t + Oi t n) ^ ] ^j,n(2) COS {p) n t T 9j r n)dt 
1 Jt ~ T n =1 n=l 

1 ft N 

— 'T' I ^ n(t) COs(u)nt T 9j, in )dj n (t') COS (uj n t 9j jn )dt 

I Jt ~ T n =1 

l ft N N 

+ — / a,, n (t) cos(oj n t + (9i >n ) ^ ay m(t) cos(uj m t + 9 jt m )dt 

n= 1 m=l, 

1 ft N i ft N 

= rfrp I ^ y Qj. njt'jO'j,nid) COS(9i, n 9j tU )dt + n(t)&?,n(t) COs(2ui n t + n + 9j tn )dt 

II Jt ~ T n= 1 ^ Jt ~ T n=l 

l ft N N 

a»,n(t)cos(a; n f+ 0*, n ) dj >m (t) cos(u; m t + 9j <rn )dt, N-A oo (3) 

n=l m—\,m^n 


Assuming that the phase of each sub¬ 
component of a chaotic carrier is a random vari¬ 
able distributed uniformly over the interval (-it, it], 
then 9i <n — 9 %n distributes over (-2-71, 2it). Hence, 
we can choose T large enough so that r %J (t) will be 
small enough. 

In CDMA systems, binary pseudo-random se¬ 
quences (chip sequences) are used to spread the 
bandwidth of the modulated signals over the 
larger transmission bandwidth, and to distinguish 
the different user signals by using the same trans¬ 
mission bandwidth. Then the chip sequence is mod¬ 
ulated by a linear sinusoidal waveform using dif¬ 
ferent modulation methods. The most commonly 
used method is QPSK modulation [Viterbi, 1995]. 
However, for simplicity and without loss of general¬ 
ity, let us choose the simplest “phase shift” method 
(“bit 1” shifts the 0 phase of the linear carrier while 
“bit 0” shifts the it phase) for demonstration pur¬ 
poses. The ith carrier in a CDMA system is given 
by 


by the kth chip value pt] namely, 


ipi(k) = 


0, p k = 1 
it, p k = 0 


Similarly, the carrier for the jth user is given 


Cj{t) — aj cos ut + 9j + E ^(*o (6) 


The moving average of the cross-correlation be¬ 
tween Ci(t ) and Cj(t) with time window T is given 

by 


l f* 

~f J T Ci ( t ) c i( t ') dt 

i r* ( [t/Tc] \ 

= — / a;cos wt+0i+ V) i>i(k) a,j 

1 J t~T \ k—O J 

/ lt/Tc 1 \ 

x cos ( ut+9j+ V'j(k) I <it 


Ci{t) = di COS Ut + 9i + ^2 


where the two constants di and to are respectively 
the amplitude and frequency of the ith carrier; T c is 
the chip bit duration defined as the time span of a 
bit in the chip sequence; [x\ denotes the biggest in¬ 
teger less than x ; and xpi(k) is the phase shift caused 


If \ 

— / didj cos 9j-9j+ y) I dt 

- 1 Jt-T / 


4/_ T a ^ COS ( : 


2 ut+9i+9j+ ' l Pi(k)+’>pj(k) dt 


( 7 ) 
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For both chaotic CDMA and linear CDMA sys¬ 
tems, the main goal is to make nj(t) —> 0 whenever 
i ^ j. Comparing Eqs. (3) and (7) we can see that 
the chaotic CDMA and the linear CDMA achieve 
this goal by using different strategies. In chaotic 
CDMA, the goal is to make the term 


I P N 

T^rr I ^ COS ( Q% n ~~ @j,n) dt (8) 

II Jt-T 

71= l 


as small as possible, but in linear CDMA the goal 
is to make the term 


1 P ( [t/Tc] \ 

2f Jt T aiajC0S I 9 i~ e o+ 1 dt 


(9) 


as small as possible. To reduce the value of 
Eq. (9) we can reduce the transmitting energy 
of each carrier (restricted by the noise level), in¬ 
crease T (restricted by the bit-rate of the mes¬ 
sage signal), increase the bit-rate of the chip 
sequence (restricted by bandwidth) and make 
tpi(k) — ipj{k) as random as possible (restricted by 
the pseudo-random algorithm for generating the 
chip sequence). 

To reduce the value of Eq. (8) we can also re¬ 
duce the transmitting energy of each carrier and 
increase T. However, instead of increasing the 
bit-rate of the chip sequence, we need to use 
as many sub-carrier components as possible; in¬ 
stead of making ipi(k) — ipj(k) random, we need to 
make 0j, n — 0 j n random enough. We can then 
conclude that the main difference between a 
linear CDMA and a chaotic CDMA is that the 
former explores the spectrum resource from the 
time-domain, while the latter does it in the 
frequency-domain. 


5. The BER Performance of (CD) 2 MA 

In this section we study the bit-error-rate (BER) 
performance of (CD) 2 MA. This is a very impor¬ 
tant benchmark for measuring service performance. 
Throughout this section we assume that all signals 
are at constant power throughout the transmission 
period. 

Suppose y % {t) is the output of the demodula¬ 
tor of the ith user and Xi >n (t) is the nth bit of 
the message signal of the ith user. As with any 
digital communication system, spread spectrum or 
not, there are four components in the demodulator 
output: 


• the desired output, which depends only on Xi >n (i)', 

• the inter-symbol interference components, which 
depend only on Xi tn+m (t), m ^ 0; 

• the component due to background noise, which 
we assume to be white with a one-side density 
equal to Nq watts/Hz; 

• the other-user interference components, which 
depend on Xj tn+m (t) for all i ^ j and all m. 

Let Ttj be the bit duration of the message signal, 
and choose a [0, T b ] integrator as the LPF. Thus, 


r = [ yi(t)dt 
Jo 


( 10 ) 


The sign of this measurement is used to decide 
whether the message bit is +1 or —1. The mean of 
T is given by 


rpb 

E[T] = [ E[yi(t)\x in (t)]dt 
Jo 

= \/Pi [ Xi >n (t)dt 


= ±T b \J~Pi 


( 11 ) 


where Pi is the power of ith chaotic carrier. The sec¬ 
ond equality is satisfied because we have assumed 
that every signal has constant power throughout the 
transmission period. Here, the sign depends on the 
sign of Xi >n (t), which is constant over [0, T b \. 

Since the noise components are essentially 
uncorrelated, 

rpb 

Var[r] = [ Yax{yi{t)\x i}n {t)]dt 
Jo 

= T b (Vi + V N + Vo) (12) 


where Vj, Vjv and Vo are the variances of the inter¬ 
symbol interference, background noise and interfer¬ 
ence from all the other users, respectively. 

Then the bit-error-rate (BER) is given by 


BER = $ 


( Km? \ 

yl Var[T] ) 





where 



(13) 


( 14 ) 
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Assuming the inter-symbol interference to be neg¬ 
ligible or zero, and since the other two components 
(background white noise and other-user interfer¬ 
ence) are independent of x^ n (t) for all i, it follows 
that the variance due to background noise is just 
the effect of a white noise with one-sided density 
No on the receiver filter whose transfer function is 
and must therefore contribute 

/ OO 

\H(f)\ 2 df = N 0 /2. (15) 

-OO 


This is because the filter gain is normalized, so that 

/ OO 

\H(f)\*df = l. (16) 

-OO 


To the user i any other user j has a constant 
carrier power P } modulated as Cj(t)xj, n (t ) which is 
independent of and generally unsynchronized with 
that of user i. In addition, the carrier phase of the 
jth user’s modulator will differ from that of the ith 
user. Hence, the effect of the j th user’s signal on 
the kth user’s demodulator will be that of a white 
noise with a two-sided density Pj passing through 
the tandem combination of two filters in the trans¬ 
mitter and the receiver with a combined transfer 
function \H(f)\ 2 . Hence, 


Vo = 



It follows from Eqs. (13), (15) and (17) that 


BER = 'f 



\ 


2 T b Pi 

\\ 

POO 

Ao + I>, / \H(f)\ 4 df 


(18) 


The numerator is just twice the bit energy, E b , and 
the denominator is the effective interference density. 
Thus, 

ber= *W?) (i9) 

where 

Eb _ _ Tt,Pi MM 

h iVo + E 

i+i 

Since (CD) 2 MA systems are interference lim¬ 
ited rather than noise limited, let us ignore the 


/ OO 

mf)\ A df 

-oo 




background noise (i.e. No = 0) and assume that 
because of perfect power control, all users are re¬ 
ceived by the base station (headend or fiber node) 
receiver at the same power P t = Pj = P c for any 
i and j. Then for a given E b /J 0 level, determined 
from Eq. (19) for the required BER, the maximum 
number of users, A max , obtained from Eq. (20), is 
given by 


Amav 1 


1 


T b 

/ OO 

\H(f)\ 4 df 

-OO 


1 

eJTo 


W/R 

OO 

\H(f)\ 4 df 

-oo 


( 21 ) 


where W is the bandwidth, R is the bit rate and 
the integral in the denominator is lower-bounded 
by unity. Consequently, 


A-may 1 


W/R 

E h /I 0 


( 22 ) 


and the maximum bit-rate, R ma x, that this channel 
can support is given by 


R 


max 


W 

E b /I 0 


+ R 


(23) 


For a given BER, the actual E b /Io depends on the 
system design and error-correction code. It may ap¬ 
proach but is never equal to the theoretical calcula¬ 
tions. In the next section, we use simulation results 
to find the performance of (CD) 2 MA systems. 


6. Comparison of BER Between 
S-CDMA and (CD) 2 MA Systems 

In the upstream, the headend receives all the 
chaotic carriers from active subscribers, additive 
impulsive noises and narrowband noises from the 
cable network. The channel model of the upstream 
data-link located at the 5-40 MHz band in cable 
TV systems is only a noisy channel with narrow- 
band interference and random impulsive noise with 
duration up to 100 p,s. Since multipath fadings 
usually encountered in wireless mobile communi¬ 
cation systems [Lee, 1989] do not occur here, a 
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n 

Fig. 8. The BER performance of the (CD) MA system used 
in a cable communication system. The dashed line and 
the dash-dotted line are the corresponding results of the S- 
CDMA system under “normal mode” (dashed line) and “fall 
back mode” (dashdotted line), respectively. 

spread spectrum communication system can be 
used much more easily in this channel than in 
a wireless communication channel. Furthermore, 
since the delay in a cable system is a fixed char¬ 
acteristic which can be measured, it is possible 
to use synchronous spread spectrum communica¬ 
tion schemes, such as S-CDMA and synchronous 
(CD) 2 MA. 

Our simulation results of the (CD) 2 MA system 
are shown in Fig. 8. The solid line shows the BER 
performance when N = 200 sub-carriers are used to 
model the chaotic carrier of each subscriber. The 
parameters for the [0, T] integrator at the receivers 
are chosen as T = T b = 1/64 ms. One should 
note that the interference from the other users in 
the same 6 MHz channel is also considered as noise 
here. For comparison, the results of the Terayon’s 
S-CDMA system in both normal mode and fall back 
mode are also shown in the same figure. We can see 
that the performance of (CD) 2 MA systems is much 
better than the S-CDMA scheme in its “fall back 
mode”, which is the best operating mode that the 
Terayon’s S-CDMA scheme can provide. 

The quantitative comparison of the channel ca¬ 
pacities of the S-CDMA in its fall back mode and 
the (CD) 2 MA is shown in Table 1. From Table 1 we 


Table 1. The relationship between the channel capacities 
of the (CD) 2 MA system and the S-CDMA system in its fall 
back mode. 

BER 10“ e 1(T 5 1<T 4 10“ 3 1(T 2 


Capacity rate: v 7 — 1.50 1.50 1.40 1.35 1.41 

CDMA 


can see that (CD) 2 MA can support a much higher 
bit-rate than S-CDMA when the BER is smaller 
than 10 -5 . A 1.5 times higher bit-rate in this case 
means that a (CD) 2 MA system can support a ca¬ 
pacity of 15 Mbps per 6 MHz upstream. 

7. Conclusions 

In this paper, we apply (CD) 2 MA to cable com¬ 
munication systems by using the 5-40 MHz band 
which is very noisy for other narrow band com¬ 
munication schemes, such as the QPSK method. 
Time division multiple access (TDMA) and fre¬ 
quency division multiple access (FDMA) also can¬ 
not be used efficiently in this noisy portion of the 
spectrum because narrowband interference and long 
duration impulsive noises can introduce many er¬ 
rors in TDMA and FDMA systems. We also com¬ 
pare our (CD) 2 MA system to Terayon’s S-CDMA 
system and found that the (CD) 2 MA system can 
perform much better than the S-CDMA system by 
increasing the channel capacity 1.5 times more than 
the best performance by an S-CDMA system. 

In this paper, we emphasize the application of 
our (CD) 2 MA on the almost “useless” spectrum 
band in today’s cable TV networks because its po¬ 
tential commercial benefit can easily be realized 
even without upgrading existing cable TV networks. 
However, readers should not form the wrong im¬ 
pression that (CD) 2 MA can only be used in coax¬ 
ial networks. In fact, the (CD) 2 MA principle can 
also be used in future high-capacity digital data 
links based on pure optical fibers to enable a high 
bandwidth efficiency. On the other hand, since 
(CD) 2 MA systems are interference limited but not 
dimension-limited, it can provide more flexible per¬ 
formance choices for different services. 
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Onset of chaotic mode transitions can be used in a simple and robust mechanism for adaptation 
in multimode systems. We present an idealized model of adaptive mode selection using on-off 
switching of chaotic basin transitions in a multistable system. A stochastic description obtained 
in the limit of large switching intervals gives search times in terms of projections of the chaotic 
state onto basins of multistability. 


1. Introduction 

This paper relates to work on harnessing the onset 
of chaotic dynamics for adaptive mode selection in 
multimode systems. This work is based on a gen¬ 
eral view of self-reorganization in which chaos can 
play an important role in finding “fit” modes. It is 
often the case that whether a particular mode of a 
multimode system, such as a laser, or a neural net¬ 
work, for example, is satisfactory or “fit”, is known 
by the external response seen when the system is in 
that mode. When working on the problem of mode 
selection in a particular high-dimensional nonlinear 
optical system, we proposed that bifurcation to in¬ 
termittent mode transitions (also called chaotic itin- 
erance) could be harnessed to automatically search 
for fit modes among a set of candidate modes by 
simply feeding back external fitness responses to a 
bifurcation parameter [Davis, 1990]. The general 
idea of adaptive mode selection using chaos is to 
couple the fitness response signal to the multimode 
in such a way that bad responses result in mode 
transitions, and good responses result in suppres¬ 
sion of mode transitions. Specifically, the fitness 
response signal can be used to drive the system 


between a multistable regime and a chaotic mode 
transition regime. 

The effectiveness and practicality of this scheme 
has been demonstrated in numerical [Davis, 1990] 
and physical experiments [Aida & Davis, 1994] us¬ 
ing an opto-electronic oscillator, and in numerical 
experiments on a neural network [Nara & Davis, 
1992]. A variation of this idea has also been used 
in experiments on a signal generator spontaneously 
adapting its output signal mode to avoid signal col¬ 
lisions [Liu & Davis, 1997]. 

The purpose of this paper is to present a simple 
model of this method, involving adaptive switch¬ 
ing between two idealized dynamical regimes. It 
allows an easy derivation of a state transition dia¬ 
gram describing the mode search process, including 
the effects of “false-alarms”. 


2. Model System 

The basic model for adaptive mode selection is as 
follows. The multimode system is described by a 
parameter p and a dynamical variable X. The 
control parameter has two particular parameter 
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values p- and p + , corresponding to the two dyami- 
cal regimes — multistability regime and mode tran¬ 
sition regime, respectively. In the multistability 
regime, there is a set of modes which are all sta¬ 
ble for the same parameter; that is, there are mul¬ 
tiple basins of attraction, and which mode is ex¬ 
cited depends on the location of the initial state. 
In the mode transition regime, the modes are un¬ 
stable, and there is a single chaotic attractor which 
extends to the neighborhoods of the previously sta¬ 
ble modes. The system receives a response E which 
is some arbitrary function of X and its history 
over some time interval. The response E takes 
two values, say 0 and 1 respectively, corresponding 
to “good” and “bad” fitness. The response signal 
causes the multimode to switch between the two dy¬ 
namical regimes at p- and p+, as shown in Fig. 1. 

The mode search algorithm is represented by 
3 steps. The fitness of the mode is evaluated, re¬ 
turning a value of signal E = 0 or 1. Then the 
parameter is adjusted, to p_ in the case of a good 
response, E = 0, or p + in the case of a bad response, 



Fig. 1. Control algorithm for adaptive mode search using 
on-off switching of chaos. It is assumed that at “test” there 
is a response signal from the environment which indicates 
whether the current state is “good” or “bad”. 


E — 1. Then the system is allowed to freely evolve 
for a time T. 

Note that we do not exclude the possibility that 
more than one mode gets a good response. And it 
is functionally reasonable for the system to make a 
spontaneous selection of any one of multiple good 
modes. 

3. State Transition Representation 

The dynamics of the search process defined in Fig. 1 
can be described by the state transition diagram 
shown in Fig. 2(a). The combined state of the mul¬ 
tistable system plus environment is completely de¬ 
scribed by the parameter p, dynamical state X, and 
response E. We sample the state of the system im¬ 
mediately after the parameter has been adjusted. 
Then there is one-to-one correspondence between 
the state of the combined system and states of X. 
First we define two states, the set Sb of X states 
which get a bad response, and the set Sc of X 
states which get a good response. One iteration 
of the search algorithm gives one of four transitions 
in the transition diagram: Sb —> Sb, Sb —> Sq, 
Sg —> Sg, or Sg —f Sb- 

Now, it is useful to further decompose Sg, as 
shown in Fig. 2(b), into “trap” states St which 
do not lead to transitions back to bad states Sb, 
and “nontrap” states Snt which do lead to transi¬ 
tions back to bad states Sb- If there is a transition 
from a good state in Sg, or any of its iterates, to 
a bad state Sb, then it belongs to the nontrap set 
Snt- The Snt states will be called “false-alarm” 
states. False-alarm states exist when there is a mis¬ 
match between the external fitness criterion and the 
basin structure. States X which get a good response 
but which are located in the basin of a bad mode 
are straightforward examples of false-alarm states. 
However, there may also be false-alarm states in 
the basins of “good” modes. Note that under cer¬ 
tain conditions, the scheme defined in Fig. 1 can 
eventually reach a good mode even if there is such 
a mismatch causing false alarms. This is a key point 
for the usefulness of this adaptation mechanism. 

4. Accessibility of Fit Modes 

Now we consider if and when the system will be 
able to find and lock onto a good mode. That 
is, do trap states exist and can they be reached 
from arbitrary initial states? If one or more of the 
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Fig. 2. (a) State transition representation of the search process. Sg and Sb are sets of states X of the multimode system 

which get “good” and “bad” responses, respectively. For this transition diagram, states are defined just after the switch 
of parameter and before the free evolution stage, (b) As for (a) but with the Sg states split into trapping states St and 
nontrapping (or false-alarm) states Snt- 


multistable modes at parameter p_ satisfy the fit¬ 
ness criterion, then the set of trap states in Fig. 2(a) 
is not empty, i.e. St # 0- If the orbit of the chaos 
at parameter p + passes through states which are in 
the trap set, then the trap set is accessible from at 
least some initial conditions. For a given system, 
such states may be easy to identify. However, in 
general, it is difficult to say whether trap sets can 
be reached from arbitrary initial states. Regardless 
of the ergodicity on the chaotic attractor at param¬ 
eter p .|_, the on-off switching may alter the distri¬ 
bution of the search dynamics to the extent that 
local accesibility is not enough to guarantee global 
accessibility. 

In order to establish a condition where we can 
guarantee global accessibility for the trap set, let 
us take the limit of large switching intervals. Let 
us assume that the interval T is much longer than 
characteristic times for relaxation to the invariant 
measures at both p_ and p + . Also, we will assume 
that each of the attractors at p- contains only ei¬ 
ther good or bad states and not both. (A counter¬ 
example is a mode which is a limit-cycle through 
good and bad states.) Then the transition diagram 
can be written as in Fig. 3, as a markov stochastic 


transition diagram. A transition between states can 
now be described as a stochastic transition charac¬ 
terized by a probability which depends only on the 
current state. The probabilities for transitions from 
the Sb state are obtained by projecting the invari¬ 
ant measure of the chaotic attractor. The probabil¬ 
ity Pb,t for transition Sb —t St is the relative mea¬ 
sure of states visited by chaos which are also in the 
trapping set St- The probability Pb,nt for tran¬ 
sition Sb -A- Snt is the relative measure of states 
visited by chaos which are in the false-alarm set. 
Now, the probability Pt,t of the St -» St tran¬ 
sition is just unity by definition. The probability 
Pnt,b °f the Snt -a Sb transition is also unity 
from the assumption about the uniqueness of the 
fitness type of the states in each multistable mode. 

It can be seen immediately from Fig. 3 that 
in this limit of large switching intervals, there is 
convergence to the trap set St so long as Pb,T is 
nonzero. It is also straightforward to calculate the 
average time r to reach the trap set, 

_ Pb,t( 1 - Pb,nt + 2 Pg ,NT + 2Pb,tPb,Nt) 

(1 - Pb,nt) 2 {Pb,t + Pb,nt ) 2 


( 1 ) 
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Fig. 3. Stochastic transition diagram obtained with assump¬ 
tion of free evolution time T longer than characteristic times 
for relaxation to invariant measures. Transitions are marked 
with corresponding transition probabilities. 


In the particular case where there are no false 
alarms, Pb,nt — 0, then the average search time 
is just 


T _ T 
Pb,T Pb,G ' 


( 2 ) 


where Pb,G is the probability of a chaotic transi¬ 
tion from a bad state to a good state, Pb,g = 
Pb,t + Pb,nt■ Clearly, the adaptive mode selec¬ 
tion process will be faster if the measure of the 
chaotic dynamics is localized on the neighborhoods 
of the candidate modes, and the environment re¬ 
sponses match the mode basins, so that the proba¬ 
bility of false alarms is reduced. Also, in practice, 
one would not want to have switch interval T any 
longer than was necessary. In this regard, for exam¬ 
ple, it would be desirable to avoid situations where 
there are good responses to states, such as those 
near fractal boundaries, which correspond to long 
transients at p_. 

By way of comparison, we could consider a 
method in which the multimode system is allowed to 
evolve in chaos at p + for time T, then uncondition¬ 
ally switched to p- where it is allowed to evolve for 
another interval T before being tested. (Note that 
in this case the response signal alone is not enough 
to switch the parameter, and supplementary param¬ 
eter dynamics are needed.) A similar analysis for 
the case of long switch intervals gives an average 


search time r of 


2 T 

T ~P^g' 

In this method, the average search time may be 
twice as long as that for the method of Fig. 1 in the 
case when Pb,nt = 0. However, this method has 
the advantage that there are cases in which it may 
allow convergence to a good mode even when the 
method in Fig. 1 does not do so. 


5. Conclusions 

We have presented a simple model of the method 
proposed in [Davis, 1990] for adaptive mode se¬ 
lection using adaptive bifurcation between multi¬ 
stability and chaotic mode transition regimes in a 
multimode system. In actual implementations of 
this method, in particular to high-dimensional mul¬ 
timode systems, such as in the experiment by Aida 
and Davis [1994], it can be difficult to analyze the 
search dynamics, for example to estimate search 
times. The simple model presented here involves 
adaptive switching between two ideal dynamical 
regimes of multistability and chaos. It allows easy 
derivation of a state transition diagram describing 
the mode search process, including the effects of 
“false-alarms”. This provides a useful framework 
for further analysis of particular systems. In par¬ 
ticular, in the limit of large switching intervals, it 
can be seen that the mode search depends on just 
the projection of the chaotic attractor on the states 
which get good responses, and on their distributions 
in the basins of stable modes. 
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A concise account is given of the early motivations for introducing parametric methods to 
achieve control of chaos. The heuristic argument that made us think that this kind of method 
could have been successful is also given. A key study is then reviewed. This concerns a 
parametric perturbation of a damped and forced Duffing-Holmes oscillator in a chaotic regime. 
The theoretical analysis, based on the Melnikov treatment of homoclinic tangles, provides a 
clear understanding of the intimate mechanism that controls chaos. Numerical results confirm 
and extend the theoretical predictions. A brief discussion of an experimental test on a magneto¬ 
elastic device is finally presented. 


1. Introduction 

It is well known that chaos is rather ubiquitous in 
both physical and nonphysical nonlinear systems. 
Sometimes chaos can be useful; this is for instance 
the case of the ergodic divertor in tokamaks, where 
a stochastic layer of magnetic field is produced at 
the plasma edge to improve the confinement. In 
other cases chaos can have harmful consequences: 
Plenty of engineering devices could be mentioned 
[Moon, 1987]. 

Among physical systems where chaos is harm¬ 
ful, we want to mention an example that motivated 
our contribution to the field of control of chaos: 
Magnetic confinement devices for controlled ther¬ 
monuclear fusion. Here the intrinsic chaoticity of 
particle dynamics is responsible for an enhanced 
diffusion across the confining magnetic field; this 
chaotic (anomalous) transport is much larger than 
the loss rate predicted by collisional transport the¬ 
ory (see e.g. [Pettini et al, 1988]). The destruction 


of regular magnetic surfaces, due to chaotic insta¬ 
bility, is another unpleasant effect in these systems 
[Rosenbluth et al. , 1966]. 

Also particle accelerators of betatron type are 
afflicted by chaotic instabilities, these can be caused 
by beam-beam interactions in storage-ring colliders 
[Scandale & Turchetti, 1991]. 

In some cases one can a priori suggest how a 
machine should be designed in order to avoid the 
onset of chaos: An example has been given for stel- 
larators [Hanson & Cary, 1984] for which the dan¬ 
gerous parameter ranges have been investigated. 

More generally, if a given physical or nonphys¬ 
ical system is satisfactorily described by some non¬ 
linear dynamical model, then by studying — ana¬ 
lytically or numerically — its parameter space, it is 
possible to know how chaos could be avoided. 

But, let us consider those situations where one 
cannot make a system operate in a safe domain of 
parameter space. In other words, assume that chaos 
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is unavoidable for the operating conditions of your 
system. For example, this is the case of anoma¬ 
lous transport of energy and particles in tokamaks. 
Then the only thing you can dream of is to perturb 
your system in a skilful way to reduce or even sup¬ 
press chaos. This idea, obviously, is not new and 
dates back to an old preprint (in Russian) [Izrailev 
& Chirikov, 1974], These authors studied how a 
perturbation of an area preserving map can change 
dramatically the phase space structure, hence the 
diffusion properties of the model; the drawback is 
in the choice of the perturbation, which is critical, 
and on its amplitude, which is not small. 

Later on, in a more recent paper [Matsumoto &; 
Tsuda, 1983], a white noise, added to a map model¬ 
ing the Belusov-Zhabotinsky reaction, was proved 
useful to reduce or suppress chaos. The explanation 
is related to the peculiar structure of the invariant 
density p(x) of the map, which is strongly peaked 
in the region of \df/dx\ that gives the largest con¬ 
tribution to the Lyapunov characteristic exponent. 
The introduction of additive noise smears out this 
peak of p(x) thus reducing chaos. 

In two subsequent papers (in Russian) 
[Loskutov, 1987; Alexeev & Loskutov, 1987], a 
parametric control of chaos has been proved effec¬ 
tive in the case of the Rossler attractor and for 
a system of ODE that models a simple ecosys¬ 
tem. These works are based only upon numerical 
simulations. 

A first account of a theoretical understand¬ 
ing of how chaos can be controlled was given in 
[Pettini, 1988], where it was presented the possibil¬ 
ity of reducing or suppressing chaos by means of 
parametric excitations on the basis of both analyti¬ 
cal and numerical results. Moreover, the suggested 
method relies upon a “resonant” effect, therefore 
a small relative variation of a parameter is effec¬ 
tive, provided that some “resonance” condition is 
satisfied. Thinking of the practical application of 
this method, its advantage is that the hardware 
of a given chaotic system should be only slightly 
modified, whereas the addition — for example — 
of new couplings in the system often requires non¬ 
trivial modifications. 

Then, at the beginning of the ’90s, an impres¬ 
sive flourishing of papers occurred on the subject of 
control of chaos [Chen & Dong, 1993]. At present 
four main strategies can be roughly identified, they 
group into two different classes: (1) adaptive — 
closed loop or feed-back — methods, requiring an on¬ 
line data acquisition and treatment of the chaotic 


signal; (2) nonadaptive — open loop or non feed¬ 
back — methods, requiring good model equations of 
the chaotic system. Both of these classes then split 
into two subclasses: (1) additive methods, where 
the orbits are directly acted upon; (2) parametric 
methods, where one or more parameters are var¬ 
ied in time in order to achieve the control of the 
dynamics. 

Claiming that one method is superior to an¬ 
other would be senseless, much depends on the 
specific problem one has to tackle. Sometimes a 
feedback method can be implemented, sometimes 
it is a priori unconceivable because following the 
individual trajectories of a system is impossible; 
this is the case, for example, of chaotic trajectories 
of charged particles in a plasma, where an open- 
loop parametric method is the only hope to control 
chaos. 

Let us now give a heuristic argument which led 
us to guess that parametric excitations could work. 
The idea arises from the following observations: 

(a) parametric perturbations can modify the sta¬ 
bility properties of fixed points of linear (or 
linearized) systems [Arnold, 1976]; 

(b) Jacobi equation for the spread of a geodesic flow 
is a linear equation whose stable and unsta¬ 
ble solutions correspond to regular and chaotic 
flows respectively. 

The first item means that the elliptic fixed point 
(x(0), x(0)) = (0, 0) of the linearized pendulum 
equation 

x T cjqX “ 0 (1) 

can be made unstable substituting lOq —» a>o(l + 
ef(t)), where f(t) = f(t + T). This is a parametri¬ 
cally excited oscillation. 

Near the hyperbolic fixed point (i(0), :r(0)) = 
(0, — n) the same equation reads 

x — loqX = 0 (2) 

and the same substitution can make stable the un¬ 
stable position (0, —7r) provided that the pivot of 
the reversed pendulum is in sufficiently rapid oscil¬ 
lation (thus e has to be large) [Arnold, 1976]. 

The second item is used only heuristically as 
follows. At least for Newtonian systems, Lagrange 
equations of motion describe the geodesics of the 
configuration space manifold equipped with the 
Riemannian metric [Pettini, 1993] < 7 ij(x) = 2 [E — 
U (x)]Jjj, where E is the total energy of the sys¬ 
tem and U(x) is the potential energy; then the 
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Jacobi equation for the second variation of the 
action functional describes the local stability of 
geodesics with respect to a reference geodesic 7 : 
{x 1 = x l (s)}; when expressed in local coordinates it 
reads [Pettini, 1993] 

Vle + R] kl x^ k x l = 0 (3) 


where Vi is the covariant derivative, Rj kl is the 
curvature tensor associated to g t j , s is the natural 
parameter along the geodesic and is the Jacobi 
field of geodesic variation. 

For two-dimensional manifolds of constant cur¬ 
vature Eq. (3) simplifies to 


ds 2 


+ A£± = 0 


(4) 


where £j_ is the perpendicular component of the 
Jacobi field £ and K is the Gaussian curvature of 
the manifold. 

From Eq. (4) it is clear that on a sphere S 2 the 
geodesics are stable because K > 0. At variance, 
on a compact hyperbolic manifold the geodesics are 
unstable because K < 0 everywhere, and thus the 
geodesic flow is chaotic. Loosely speaking, to de¬ 
scribe regular and chaotic dynamics we have recov¬ 
ered — at another level — Eqs. (1) and (2). Let¬ 
ting K —» K(1 + ef (t)), as with Eq. (1), one can 
make exponentially unstable nearby geodesics on a 
positively curved manifold as a consequence of cur¬ 
vature fluctuations “felt” by the geodesics; this is 
actually a major mechanism responsible of chaos 
in Hamiltonian flows of physical relevance [Cerruti- 
Sola Sz Pettini, 1996; Pettini & Valdettaro, 1995; 
Casetti et al, 1996]. Therefore it is also conceivable 
that a suitable parametric perturbation of Eq. (4) 
might act to stabilize the exponentially unstable 
(chaotic) trajectories, when K < 0, in analogy with 
the reversed pendulum Eq. (2). Within this analogy 
the sign of K should periodically change in time. 

This argument is only heuristic because, in gen¬ 
eral chaotic flows are not topologically equivalent 
to geodesic flows on manifolds of constant nega¬ 
tive curvature, if this were the case one should 
have structural stability (after the Lobatchevsky- 
Hadamard theorem [Arnold & Avez, 1968]), thus 
ergodicity, mixing, etc., but this is not the generic 
situation. In principle this heuristic argument — 
of differential geometric kind — could be also re¬ 
peated for dissipative systems using the geometry 
of Finsler spaces, but this goes far beyond the aim 
of the present contribution. 


The conjecture given above is tested using a dy¬ 
namical system where, to some extent, chaos can be 
tackled also analytically. 


2. A Paradigmatic System 

In [Lima & Pettini, 1990] we chose the so-called 
Duffing-Holmes oscillator. This model, defined by 
the equation 


x — ax + (3x 3 = — Sx + 7 cos ut (5) 

is one of the simplest nonlinear dissipative ODE un¬ 
dergoing a chaotic transition. With some approxi¬ 
mations of Galerkin type [Guckenheimer & Holmes, 
1983], it can be derived from a PDE describing the 
dynamics of a buckled beam; in a different context, 
it can also be used to describe plasma oscillations 
[Laval & Gresillon, 1979]. Equation (5) can be triv¬ 
ially rewritten as (we set a = 1 without loss of 
generality) 

(*) = ( V B 3 )+e( 5 7 ) (6) 

\y) \x - fix 3 ) + ^cos ut) 

which is in the form 


x = f 0 (x) +efi(x, t). (7) 


The unperturbed part x = fo(x) can be derived 
from the Hamiltonian 

H = \ y2 ~ \ x<1 + ( 8 ) 

and is integrable. Its phase space has only one hy¬ 
perbolic fixed point from which an “eight-shaped” 
separatrix originates. The motion on this separatrix 
is given by 


a S°\t) = y — secht, 

p2 

x(°)(t) = —J—secht tanh t. 


(9) 


The separatrix, parametrically defined by Eq. (9), 
is also called homoclinic loop and results from the 
superposition of the so-called stable and unstable 
manifolds, W s and W u , respectively tangent at the 
origin to the stable and unstable eigenspaces E s and 
E u of the hyperbolic point. W s and W u are defined 
as those trajectories which converge asymptotically 
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to the hyperbolic fixed point: W s for £ —> +oo 
and W u for £ —»■ —oo respectively. When the sys¬ 
tem x = fo(x) is perturbed only by a dissipative 
term, the two manifolds W s,u never meet and the 
solutions are still regular. If a forcing term is also 
added (i.e. an energy supply is added to balance 
friction losses) then W s,u may have an homoclinic 
intersection and hence an infinity of subsequent in¬ 
tersections [Guckenheimer & Holmes, 1983]. We 
briefly recall how Melnikov’s method works to de¬ 
termine the condition of homoclinic intersection of 
W s and W u and so of the onset of chaos. Let 
rW(t) = (a;( 0 )(£), x(°\t)) T be the unperturbed mo¬ 
tion on the homoclinic loop, write 

w s ’ u {t, to) ~ r <°>(£ - to) + sW s,u i(t, to) (10) 

to describe how W s ' u are perturbed up to first order 
in e [due to fi in Eq. (7)] starting from r(°); to is an 
arbitrary reference time and W s,u = ( x s,u , x s ’ u ) T 
are column vectors. One gets 

= J(r (0) (t - t 0 ))IE s - u i 
+ ( 11 ) 
where J is the Jacobian matrix of fo computed at 

r<°> (£-£«>). 

Then the Melnikov distance is defined as 

A (£, to) = n • to) - W?(t, to)) (12) 

where n is the normal to T^(t — to)- 

After some algebra one finally finds the Mel¬ 
nikov function 


/ OO 

dt(fo A fi)p(o)( t _ to) (13) 

-OO 

which in principle can be explicitly computed; if 
A (to) changes sign for some to, then an infinity of 
homoclinic intersections between W u and W s will 
take place and chaos will set in. This is the only 
general predictive method to study the condition 
for the onset of chaos in dissipative ODE. 

Notice that for Hamiltonian systems there are 
always homoclinic intersections when a noninte- 
grable perturbation efi(x, t) is added to an inte¬ 
grate system; in this case the Melnikov function 
[Chirikov, 1979] is 


/ OO 

dt{Ho, #i}r(o)(t-to) i (14) 


where curly brackets are Poisson brackets, of the 
unperturbed Hamiltonian Hq with the perturba¬ 
tion Hamiltonian Hi, computed along the unper¬ 
turbed separatrix I^ 0 ); M(to) is useful to evaluate 
the thickness of the stochastic layer. The analytical 
computation of A(£o) for Eq. ( 6 ) is standard and 
yields 


A(to) = 27 ta —qcjsech 


TTUJ \ 4 8 

Tj sinwi ° + 5r (is) 


Unfortunately there are not so many models for 
which explicit computation of A(£ 0 ) can be per¬ 
formed. Therefore we chose Duffing-Holmes model 
because it is not difficult to compute A(io) when a 
parametric perturbation is introduced. 

Let us modify Eq. ( 6 ) to 


_ / v 

yj yx — (3( 1 + r) cos £lt)x 3 


+ £ 


5 7 

—V H— cos cot 
£ £ 


(16) 


if 77 <C 1 we are allowing a periodic modulation of 
small amplitude of the parameter (3. 

Accounting for the modulating term 
/??7 cos fit x 3 in the perturbation function fi, the 
new Melnikov function A p (£q) is given by 


A p (t 0 ) = A (£ 0 ) - -If dtx^(t - t 0 ) 

P J — OO 

x[x(°)(£-£ 0 )] 3 cos fit (17) 


using (10), after simple but tedious computations, 
one finds [Lima 8z Pettini, 1990, 1993] 


. . . „/2 , / 7ru;\ . 46 

A p (£ 0 ) = 2 ttW - 7 u> sech j sin u >£ 0 + ^ 

— w^(D 4 + 4D 2 )cosech f —^ sin flto - 
up \ 2 / 

(18) 


Let us now consider a set of parameters for which 
A(£o) changes sign, thus predicting the existence 
of chaos; then we add the parametric perturbation 
as in Eq. (16) and we compute, from the numerical 
tabulation of Eq. (18), the time needed between two 
successive homoclinic intersections of W u and IV s , 
let us denote it by tm- 
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Fig. 1 . The inverse of the time tm elapsed between two suc¬ 
cessive homoclinic intersections, computed using Eq. (18), 
is plotted versus the parametric perturbation frequency fi. 
Again (3 = 4, 5 = 0.154, 7 = 0.088, oj = 1 . 1 . Continuous line 
corresponds to 77 = 0.09, dashed line to r] = 0.1 and dotted 
line to r) = 0.15. 


In Fig. 1, r ^ 1 is reported versus the paramet¬ 
ric perturbation frequency XI for the following set of 
parameters: 7 = 0.088, S = 0.154, u — 1.1, /? = 4, 
77 = 0.1. 

A “resonance line” is found which is centered 
at XI = oj] for this value of XI, tm becomes infinite; 
this means that W u and W s never intersect, hence 
chaos should disappear. 

If 77 is smaller than some critical value r] c , then 
tm remains finite and homoclinic intersections are 
not suppressed. By increasing 77 above T] c , a line 
broadening is observed (see Fig. 1). Finally, outside 
the interval of SI where the quartic polynomial in 
Eq. (18) is negative, there is no way to avoid inter¬ 
sections of W u and W s . 

Similar results are obtained for different sets of 
parameters. 

Accurate numerical experiments can be per¬ 
formed to make a comparison with these pre¬ 
dictions. Equation (16) is integrated using a 
Hammings modified predictor-corrector of fourth 
order, time integration steps At = 0.001-0.003 and 
integration times of t = 10 000 — 20 000 after a tran¬ 
sient of t = 500. By means of a standard technique 
[Benettin et al. , 1979], the largest Lyapunov char¬ 
acteristic exponent A is computed to detect chaos 
and measure its strength. 



Fig. 2. The largest Lyapunov characteristic exponent A 
is reported versus Q. This result corresponds to ft = 4, 
S = 0.154, 7 = 0.088, oj — 1.1, rj = 0.03. The dotted line 
refers to 77 = 0 (unperturbed case) for which A = 0.1056. 
Dotted lines are also used for subharmonic resonances whose 
existence depends on the starting point in phase space. 

For the same set of parameters in Fig. 1, ex¬ 
cept for 77 which is 0.03, A is computed for differ¬ 
ent values of the parametric perturbation frequency 
H. The results are shown in Fig. 2 and are rather 
striking. The resonance line displayed in Fig. 1 by 
t m (^)> ^ Xl° R = 1.1 corresponds to the first order 
resonance of X(Xl) in Fig. 2. Moreover, other reso¬ 
nances show up: at XI — 2.2 and XI = 3.3, second 
and third harmonics of Xl° R respectively, at XI = 0.55 
and 12 = 2.75, i.e. at the first subharmonic of 
Xl^ and at the third harmonic of the subharmonic. 
These last resonances are very narrow and corre¬ 
spond to suppression of chaos (A = 0) only if the 
initial conditions belong to a suitable domain. 

In Fig. 3 it is shown that a resonance broad¬ 
ening is produced at increasing perturbation ampli¬ 
tude, which is at least in qualitative agreement with 
the analytical result reported in Fig. 1. 

Another way to get a hold of what happens 
to the dynamics when the parametric perturbation 
frequency approaches a resonant value, is to look at 
the autocorrelation function (x(t + A t)x(t)) of the 
solution of Eq. (5). In [Lima & Pettini, 1990] some 
of them are reported and show an increasing corre¬ 
lation time — thus a decreasing chaoticity — when 
Xl approaches a resonant value (in this case X^' R ). 

We thus have a paradigmatic example show¬ 
ing that chaos can be reduced or eliminated in a 
dissipative system by means of parametric pertur¬ 
bations. A 3% modulation of a parameter, when 
the modulation frequency is resonant with the forc¬ 
ing frequency, is able to make regular the chaotic 
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Fig. 3. A versus Cl for /3 = 4, S = 0.154, 7 = 0.088, uj — 1.1 
at different parametric perturbation amplitude. Dotted line 
corresponds to r) = 0.03, continuous line to r) = 0.05 and 
dashed line to rj = 0.07. 

dynamics. Let us just spend few words about 
Hamiltonian systems. We have already mentioned 
that in this case chaos is always present, therefore 
one can “only” hope either to stabilize some bro¬ 
ken KAM tori, or to modify the diffusion proper¬ 
ties of the model. Some attempts made with the 
Hamiltonian version of the Duffing-Holmes oscilla¬ 
tor and with the systems of [Pettini et al., 1988], 
show that sometimes several small and sticky is¬ 
lands can be produced by parametric perturbations; 
the major consequence is a slowing down of diffusion 
in phase space due to intermittent trapping near is¬ 
lands. An increase of intermittency is revealed by a 
worse convergence of maximal Lyapunov exponent 
caused by enhanced fluctuations of local divergence 
rate of nearby trajectories. In general, stronger 
perturbations are necessary to produce measurable 
effects. 

3. An Experimental Confirmation 

The model equation above studied is well suited also 
for an experimental test of the practical applicabil¬ 
ity of the method [Fronzoni et al., 1991]. 

The experimental apparatus consists of a 
magnetoelastic device schematically represented in 
Fig. 4: A steel-made elastic beam is clamped at 
one end and is left free at the other end; near the 
free end of the beam two magnets create a two- 
well potential with two stable equilibria. The whole 
system is put in vibration by means of an electro¬ 
magnetic shaker. If the forcing is not too strong, 



Fig. 4. Sketch of the experimental apparatus. 


a one mode approximation of Galerkin type can be 
done on the partial differential equation that de¬ 
scribes this flexible pendulum [Marsden & Holmes, 
1980], so that the single mode amplitude obeys the 
Duffing-Holmes ordinary differential Eq. (5). In or¬ 
der for a parametric perturbation to act upon the 
system, the permanent magnets of the experimental 
apparatus are surrounded by coils of copper wires. 
An oscillating electric current in the coils gives a 
modulation of the magnetic field at the end of the 
beam so that a and (3 in Eq. (5) become 

a —¥ a[l + e cos(27n^f)], £<1 (19) 

/?-»/?[I + 77 cos{2m/ M t )], 77 < 1. (20) 

The bending of the beam can be measured by 
an optical device: A thin screen fixed on the elastic 
beam makes an occultation of a light source, then 
the resulting light intensity, measured by a photore- 
sistence, reveals the fundamental mode amplitude 
of the deflection which is described in Eq. (5). 

The dynamics of the system can be examined 
with the aid of an oscilloscope where the portraits 
{x, x ) of the Poincare sections of the phase space 
are obtained by modulating the z-axis of the oscil¬ 
loscope with a signal synchronized to the driving 
signal of the vibrator. A typical Poincare section, 
obtained when chaotic vibrations are present, is 
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Fig. 5. Poincare section obtained with 300 seconds of acqui¬ 
sition time. About 1000 points are displayed, uf = 13.03 Hz, 
A/A c = 1.33, e = Tj = 0. 


reproduced in Fig. 5. Long runs (several hours) 
of the experimental apparatus ensure that chaos is 
stable in the system. 

A current flowing through the coils around the 
permanent magnets changes the magnetic field in¬ 
tensity at the end of the elastic beam, this means 
that the parameters a and f3 are modified according 
to Eq. (20). 

When a modulation voltage is applied to the 
coils, the response is like in Fig. 6, where the volt¬ 
age output (proportional to x(t)) is displayed ver¬ 
sus time. At modulation frequencies I'm close to 
the forcing frequencies up chaotic vibration is al¬ 
ternated by ordered vibrations. Ordered phases are 
almost periodic oscillations around one of the two 
minima of the potential with x(t ) locked to the forc¬ 
ing signal. 

The average duration of the regular vibrations 
increases when \vp— vm\ is reduced. Roughly speak¬ 
ing, the parametric perturbation makes the system 
intermittent and by approaching a resonance be¬ 
tween up and vm the laminar phase has an increas¬ 
ing weight with respect to the chaotic phase. Com¬ 
plete regularization of the dynamics shows up only 
at exact resonance. A modulation of the 10% of the 
unperturbed value of a (i.e. e = 0.1) is needed to 
produce the result reported in Fig. 6. 

A less qualitative description of phenomenon is 
given in Fig. 7, where the amplitude output of a 
Lock-in measurement device is reported versus the 
perturbation frequency I'm when the forcing signal 



Fig. 6. Oscillation amplitude x(t), measured in Volts, at 
u F = 13.03 Hz, e = 0.1, 7? = 0.08, A/A c = 1.33. 
(a) vm — 13.06 Hz; (b) um = 13.07 Hz; (c) vm = 13.08 Hz. 


is sent also to the reference of the Lock-in. The 
tracings display typical resonance patterns. These 
show how the relative weight of laminar phases 
with respect to chaotic phases changes as a func¬ 
tion of parametric modulation frequency. The two 
upper curves correspond to different amplitudes 
of the modulation (e = 0.16 and e = 0.1 from 
top to bottom) and show a resonance broadening 
at increasing e. The maxima of the peaks corre¬ 
spond to I'm — vp. The same effect of parametric 



1682 R. Lima M. Pettini 



28 30 V m (Hz) 


Fig. 7. Lock-in amplitude output (arbitrary units) versus 
I'm- Sweep time: 1175 sec/Hz. Lock-in integration time: 
300 sec. vp — 14.59 Hz and A/A c = 1.33. (a) upper curve 
refers to e — 0.16; lower curve refers to e = 0.1. (b) second 
harmonic resonance, e = 0.12. 


modulation is found at vm — 2 uf and this too is 
reported in Fig. 7 (lower curve). 

In the experimental device it is not possible to 
avoid the contemporary modulation of a and 0, 
therefore one has to adapt the above theoretical 
analysis to this case. The Melnikov function must 
be computed for Eq. (5) with the modulations given 
by Eq. (20). The new result — (setting $2 = 2TruM 
and uj = 2m/p) — reads 


A (to) = 2nu>A 




7rea 2 

(0! d 

. 0 

U ) 


7TT] 

6/3 


(fl 4 +4fl 2 ) 


x cosech ( sin(fiio) + 

,2^/aJ 60 


( 21 ) 


At 'q = £ = 0, if the parameters are such that 
the function A(/o) changes sign then homoclinic in¬ 
tersections are present and hence the existence of 
chaotic solutions is inferred. At £,rj ^ 0, as the 
sign of the polynomial function of S2 is indefinite, 
the second term of Eq. (21) can act as a counterterm 
of the first one. As a consequence, also choosing the 


unperturbed system in a chaotic phase, it is again 
possible to prevent A(fo) from changing sign (i.e. to 
eliminate chaos), provided that Q is resonant and 
that the other parameters are in suitable domains. 


4. Concluding Remarks 

We have briefly shown that open-loop parametric 
control of chaos is successful. The model here dis¬ 
cussed is somehow paradigmatic because, beside the 
general heuristic argument, we have at disposal an 
analytical prediction — based on the only predic¬ 
tive method existing at present (Melnikov’s theory) 
— and numerical simulations confirming the theo¬ 
retical analysis, and finally an experimental confir¬ 
mation on a physical system. 

Numerical simulations confirmed the analyti¬ 
cal predictions revealing at the same time a richer 
phenomenology, which is not surprising because 
Melnikov’s theory is based on an approximate 
method. 

The experimental confirmation is crucial at 
least for two reasons: The obvious one is to show 
that the method can be successfully implemented 
in real systems, the other is that the possibility of 
controlling chaos by parametric perturbations turns 
out to be robust, in fact our experimental system is 
only roughly modeled by a Duffing-Holmes equa¬ 
tion, nevertheless the method works, this suggests 
that the choice of good model equations is certainly 
important but not critical. 

Many problems are open and need a deeper un¬ 
derstanding. Among the open questions we mention 
that: 


(i) The development of a more general theoretical 
framework requires to go beyond Melnikov’s 
theory which is undoubtely general but it ap¬ 
plies to weakly perturbed systems and it is of 
practical utility only in a few simple cases. An 
interesting direction for further research work 
on this topic could be provided by optimal con¬ 
trol theory, a la Pontrjaguin, i.e. to minimize 
the functional 


J = 



dtA(x(t), u(t )) 


rt\ _j__ 

+ 0 dt ^ - /i(x(t), U(t))} 2 (22) 

Jt ° i =i 

where the dynamic constraints iq = fi(x(t), 
u(t)) have been introduced with a Lagrange 
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multiplier p, and where the (vector) function 
u(t) represents the control(s) with respect to 
which the functional J has to be minimized. 
Here the problem is to find a good function 
A(x(t), u(t)) that attains its minimum in cor¬ 
respondence of ordered motions. 

(ii) The above results raise an intriguing prob¬ 
lem in connection with the so-called “stabil¬ 
ity dogma” [Guckenheimer &; Holmes, 1983]: 
If you observe chaos in real systems it must be 
structurally stable (in the strong sense) and 
the same should hold for the theoretical mod¬ 
els used. In fact theoretical models will never 
take into account all the interactions, pertur¬ 
bations, noise, etc. which are present in real 
systems. Attempts to soften the definition of 
structural stability have been proposed some 
years ago [Zeeman, 1988]. 

(iii) In this context it is not out of place to men¬ 
tion that, after the Birkhoff-Smale homoclinic 
theorem [Guckenheimer & Holmes, 1983], the 
existence of homoclinic intersections for the 
Dufhng-Holmes oscillator ensures the exis¬ 
tence of a hyperbolic invariant set A. Other 
considerations [Guckenheimer & Holmes, 1983] 
rule out, for the same model, the possibility for 
A to be an attractor. There are several reasons 
to believe that noise, which is always present 
both in real systems and in numerical mod¬ 
els, plays an important role together with A 
to stabilize chaotic transients or, at least, to 
make them very long with respect to practical 
observational times. 
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This paper describes mode-transition dynamics of a compound-cavity, multimode, Fabry-Perot 
(FP) semiconductor laser and the basic idea of adaptive wavelength selection using chaotic 
search. In the presence of external feedback, the laser is shown to exhibit different mode- 
competition dynamics including single-mode, multistability, and chaotic mode-transition states. 
A chaotic search technique is proposed to adaptively select a dominant lasing mode which fits 
the environment by switching a control parameter (feedback level or bias injection) between 
multistability and chaotic mode-transition states. The effectiveness of the adaptive mode selec¬ 
tion is verified through a computer simulation and the robustness as well as the experimental 
implementation of the technique are discussed. 


1. Introduction 

Applications of chaotic dynamics have received 
much attention in recent studies of nonlinear dy¬ 
namics. One promising topic is adaptive mode se¬ 
lection using chaotic search [Davis, 1990] which uti¬ 
lizes the global properties of a chaotic system and 
the bifurcation from multistable regime to chaotic 
regime. This method was successfully applied to 
adaptive oscillation mode selection in an electro- 
optical system [Aida & Davis, 1994; Liu & Davis, 
1997], 

This paper describes the application of the 
chaotic search to a multimode FP semiconductor 
laser with external optical feedback. Experiments 
have revealed that the external reflectivity as well 
as the external cavity length exert essential influ¬ 
ences on dynamical and spectral behaviors of a laser 
diode in an external cavity (e.g. [Mprk et al., 1992; 
Uenishi et al., 1996]). The purpose of this paper 
is to investigate the mode-competition dynamics 
in the external cavity laser and to verify whether 
it has a bifurcation structure which is suitable 


for performing adaptive mode (wavelength) selec¬ 
tion using chaotic search, that is, multistability of 
modes and an onset of chaotic transitions among 
the different modes. 

The paper is organized as follows. In the next 
section, we first introduce a numerical model of a 
multimode laser with external feedback and com¬ 
pare numerical results for the spectral distribu¬ 
tion with experimental observations. The mode- 
competition dynamics are characterized by mode 
power and dominant mode distributions, and based 
on these measures, the laser output is classified 
into three different regimes: single dominant mode 
regime (only one mode dominates the lasing), multi¬ 
stability regime (different modes become dominant 
depending on initial condition), and chaotic mode- 
transition regime. In particular, we show that it 
is possible to switch between the multistability and 
chaotic mode transition regimes by switching a con¬ 
trol parameter: external feedback strength or injec¬ 
tion current. In Sec. 3, we show that adaptive mode 
selection using chaotic search can be performed by 
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coupling the control parameter to an external eval¬ 
uation signal. The effectiveness of the adaptive 
mode selection technique is verified by computer 
simulations. The robustness and possible exper¬ 
imental implementation of the technique are also 
discussed. The last section summarizes the paper 
with a comment on potential applications of the 
proposed chaotic search scheme to other lasers. 


2. Mode-Competition Dynamics and 
Bifurcation Diagrams 

As a reference for the numerical modeling, we first 
show an experimental optical spectrum of a typ¬ 
ical laser diode with moderate optical feedback 
(1 ~ 10%). The laser employed in the observa¬ 
tion of Fig. 1 is a AlGaAs laser diode with a center 
mode wavelength of about 782 nm at the threshold 
bias current of 45 mA. It is well known that a FP 
semiconductor laser subject to moderate to strong 
feedback often exhibits multimode behavior. In this 
experiment, it is found that about 10 modes have 
power levels within 10 % of that of the peak mode 
power. 

The dynamics of the above laser can be de¬ 
scribed by a set of multimode rate equations in¬ 
cluding the external feedback and Langevin noise 
[Byrne, 1992; Ryan et al., 1994]. 


dE m (t) .. , . 

dt = - Wth )Em{t) 

+ 2 (! “ ia )(Gm - 7 )E m (t) 

+ KE m (t - t) exp (iu) m T) 


- I £ 9 mk \E k (t)\ 2 E m (t) + F m (t ), ( 1 ) 

Z k 


dNjt) 

dt 


m 

T s 


M 

-J2 G j\ E M 2 + F N (t). 


j =i 


( 2 ) 


Here, E m is the complex optical field of the mth 
(m = 1,2,..., M) longitudinal mode and is scaled 
so that \E m \ 2 corresponds to the photon number 
within the laser cavity, N(t) is the carrier number 
within the laser cavity, a> m , A m , G m are respectively 
the angular frequency, wavelength, and gain of the 
mth mode. 7 is the loss coefficient of the laser, r s is 
the life time of carriers, r is the time delay in the ex¬ 
ternal feedback, oj t ,h is the angular frequency of the 
solitary laser without external feedback, a is the 



Wavelength (nm) 
(a) 



Mode 


(b) 

Fig. 1. Optical spectra of a FP laser diode with external 
feedback, (a) Experimental and (b) simulation results. The 
vertical scale of (a) is 50 nW/div. 


linewidth enhancement factor, 0 m fc is the coefficient 
for cross-saturation between modes m and k, J is 
the bias injection, and F m , Fn are Langevin noise 
terms. We did not include noises in the current sim¬ 
ulations in order to distinguish effects of determin¬ 
istic chaos from stochastic noise. The expression for 
the mode gain G m is given by 

G m = g[N{t) - Nq \[1 - (A m - A p ) 2 /AA;;], (3) 

where g is the gain coefficient, No is the carrier num¬ 
ber at transparency, A p is the wavelength of the 
peak gain, and AX g is the full width half maximum 
(FWHM) value of the parabolic gain curve. The 
feedback level k is given by 

k = (1 - ro)ri/r in r 0 , (4) 

where Tj n is the round-trip time of the laser cav¬ 
ity and ro and r\ are internal cavity reflection 
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Table 1. Some parameter values for the laser diode used in the 
numerical calculations. 


Symbol 

Parameter 

Value 

A 

Center wavelength 

780 nm 

SXm 

Wavelength separation 

0.28 nm 

SXg 

FWHM of laser gain 

100 nm 

1 

Laser cavity length 

250 mm 

ro 

Facet reflectivity 

0.5556 

9 

Gain coefficient 

1639 s _1 

7 

Loss coefficient 

5.3 x 10 11 s _1 

a 

Linewidth enhancement factor 

3 

No 

Carrier number at transparency 

1.75 x 10 8 

N t h 

Carrier number at threshold 

4.99 x 10 8 

T s 

Life time of carrier 

2.0 ns 

7*in 

Round-trip time of laser cavity 

5.9 ps 

@mm 

Saturation coefficient 

2.8 X 10 4 s _1 

@mk 

Cross-saturation coefficient 

2.5 X 10 4 s -xl 


and external reflection coefficients (both in ampli¬ 
tude), respectively. Modes couple with each other 
through the depletion of the commonly shared car¬ 
riers [Eq. (2)] and the gain cross-saturation effect. 
We neglect the small contribution from spontaneous 
emission to the lasing mode. This factor is evalu¬ 
ated to be much smaller (less than one percent) than 
the feedback term (at ri = 1%) and negligible com¬ 
pared with both the gain and loss factors. Based 
on the observed optical spectrum, we assume the 
bandwidth of the laser gain to be 100 nm with the 
center wavelength at 780 nm and the total mode 
number M to be 20 which is large enough to cover 
the experimental spectrum. Other parameters and 
parameter values employed in the calculation are 
listed in Table 1 [Lang & Kobayashi, 1980; Mprk 
et al., 1992; Ryan et al ., 1994]. 

We have numerically simulated Eqs. (1) and 
(2) by employing a fourth order Runge-Kutta al¬ 
gorithm and investigated dynamical states for dif¬ 
ferent parameter conditions. Figure 1(b) shows a 
numerical result for the optical spectrum of the laser 
with external feedback. The spectrum is quite sim¬ 
ilar to the experimental result in Fig. 1(a) and this 
justifies the numerical model and the parameter val¬ 
ues used. 

To characterize mode states of the multi- 
mode LD, we introduced two measures, namely, 
mode-power-ratio (MPR) and dominant-mode- 
ratio (DMR). MPR stands for a long-term average 


of the light power of each mode. The distribution of 
MPR shows the optical spectrum of the laser during 
a specified time interval. Meanwhile, DMR shows 
how frequently a specific mode dominates the laser 
output. The two factors are respectively defined as 

(MPR)j = j lE^ttfdt/Y, J \Ej(t)\ 2 dt, (5) 

and 

(DMR)j = /|n W)-^(«)]W/ dt, (6) 

l j^i ) ' 

where 9{x) is a step function and takes 1 for x > 0 
and 0 otherwise. 

We mention that nonlinear dynamics and mode 
interactions of multimode lasers have been stud¬ 
ied in various experimental systems and numer¬ 
ical models. For example, Byrne [1992] calcu¬ 
lated the averaged response and the stochastic 
response of a multimode laser under Gb/s mod¬ 
ulation. Bracikowski and Roy [1991] investigated 
periodic pulsing of a multimode solid-state laser 
and found antiphase cycling among different modes. 
Ryan et al. [1994] investigated intensity noise fea¬ 
tures of a multimode laser diode with external feed¬ 
back and observed mode transitions in the time 
evolution of the laser output. Recently, Szwaj 
et al. [1996] observed wave propagation in the 
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Single Mode Multi-Stability Mode-Transition 

Regime Regime Regime 



External Reflectivity (%) 


Fig. 2. Bifurcation diagram of mode-competition states ver¬ 
sus external reflectivity. White circles and solid dots denote 
respectively stable and temporary dominant modes. Param¬ 
eters are h = 1.5 7 t h and L ext = 15 cm. 


klode-Transition Multi-Stability Mode-Transition 

i Regime Regime Regime 



Bias Injection (Ith) 


Fig. 3. Bifurcation diagram of mode-competition states ver¬ 
sus bias injection. White circles and solid dots denote respec¬ 
tively stable and temporary dominant modes. Parameters are 
Lext = 36 cm and R ex t = 1.0%. 


spectrum as a result of coupling between modes in 
a strongly multimode fiber laser. In this paper, we 
focus on the bifurcation of mode-competition dy¬ 
namics when the external feedback level or the bias 
injection are varied. In particular, we use the mea¬ 
sures of DMR and MPR to classify different laser 
states and investigate how the bifurcation of mode- 
competition dynamics results in changes in laser 
states. 

Figures 2 and 3 show the bifurcation diagrams 
of laser states versus the external reflectivity and 
the bias injection, respectively. In these diagrams, 


white circles denote stable dominant modes whose 
DMR is unity for some initial condition. Multi¬ 
ple circles at a fixed parameter value indicate mul¬ 
tistability of the dominant mode. On the other 
hand, solid dots denote temporary dominant modes, 
i.e. the dominant mode changes from one to another 
in time and each marked mode is only dominant for 
a short time interval. It is stressed that the blank 
area in the bifurcation diagram does not necessarily 
mean the output power of the mode there is truly 
zero, it only means the mode never becomes domi¬ 
nant during the observation time. 

In the above bifurcation diagrams, the pa¬ 
rameter domain can be divided into three dif¬ 
ferent regimes corresponding to different mode- 
competition states: (i) single mode regime where 
only one specific mode becomes dominant (which 
means DMR = 1) and the mode takes most of the 
output power all the time; (ii) the multistability 
regime where different dominant modes exist de¬ 
pending on initial conditions; and (iii) the chaotic 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
Mode 

(b) 

Fig. 4. Mode distributions in multistability regime for R = 
1.5 7 t h, Lex t = 15 cm, and 7? ex t = 1.4%. Dominant mode is 
(a) No. 9 and (b) No. 13. 
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mode-transition regime where no particular mode 
becomes dominant during the observation time. 

Examples of mode distributions for multistable 
states are shown in Figs. 4(a) and 4(b). The his¬ 
tograms show DMR distributions while the lines in¬ 
dicate MPR values. In both figures, there exists 
one mode whose DMR becomes unity, although the 
MPR shows that the mode does not take all of the 
laser output power. 


Thus the search is performed as follows. When a 
mode gets a good external response, the control pa¬ 
rameter is set to /iQ. If the laser output falls into a 
different mode which gets a bad response, the con¬ 
trol parameter is reset to /z c and the search loop 
is repeated. 

Two typical examples of chaotic search are 
shown in Figs. 5 and 6. In the case of Fig. 5, the 


3. Adaptive Mode Selection Using 
Chaotic Search 

A method for adaptive mode selection using chaotic 
search was proposed by Davis [1990]. Here the word 
“mode” could mean an oscillation pattern, a bit se¬ 
quence, or the dominant lasing wavelength as in the 
current case. The method is applicable to systems 
which exhibit multistability of modes and an onset 
of chaotic transitions among the different modes. 
It is further assumed that the “fitness” of a partic¬ 
ular mode is shown by an external response signal 
obtained when the system is in that mode, with dif¬ 
ferent response signal levels corresponding to good 
fitness and bad fitness. The general idea of adap¬ 
tive mode selection using chaotic search is to couple 
the fitness response signal to the control parameter 
of the system in such a way that bad responses re¬ 
sult in mode transitions and good responses result 
in mode stability. 

As shown in the last section, a multimode FP 
laser diode with external feedback has multistability 
and chaotic mode-transition states as well as bifur¬ 
cation between these two regimes upon the contin¬ 
uous variation of laser parameters. Thus, we are 
able to apply the chaotic search idea to the mul¬ 
timode FP semiconductor laser for the purpose of 
adaptive mode selection. The chaotic search con¬ 
sists of two steps: mode evaluation and parameter 
adjustment. In the first step the system output is 
evaluated and external response is generated and in 
the second step the control parameter is adjusted 
using the external response signal. In the first step, 
we detect the output of each mode and calculate the 
DMR during the fixed mode-evaluation time inter¬ 
val. The mode-evaluation time is chosen as several 
times the relaxation time of the laser. For the sec¬ 
ond step, the control parameter is chosen as a two- 
level variable for simplicity, with parameter value 
/i c corresponding to chaotic mode-transition regime 
and parameter value p.° to multistability regime. 


13 


13 


li 



Fig. 5. Mode search process. Target modes are mode 8, 
13, 11, and 13. The control parameter (bias injection /&) is 
switched between p c = 1.15 7th and p° = 1.35 7 t h- Other 
parameters are L ext = 36 cm and 7? e xt = 1.0%. Upper trace: 
Dominant mode DMR. Lower trace: The control signal. 



Time (|is) 


Fig. 6. Mode search process. Target modes are mode 13, 9, 
11, 9, and 13. The control parameter (external reflectivity 
/text) is switched between p c = 2.8% and p° = 1.4%. Other 
parameters are A = 1.5 7 t h and L ex t = 15 cm. Upper trace: 
Dominant mode DMR. Lower trace: The control signal. 
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bias injection is taken as the control parameter 
and is switched between chaotic mode-transition 
regime at /x c = 1.15 7 t h and multistability regime at 
/x° = 1.35 7 t h- Three modes (mode 8, 11, and 13) 
were alternatively set as targets. During each search 
loop, there are several switch-on and -off processes. 
Such behavior is typical in this adaptive mode selec¬ 
tion method due to mismatch between the external 
classification of output states and the basins of at¬ 
traction. When the target mode appears, the laser 
is set in multistability regime (p = /r°) and the out¬ 
put mode converges to one of the multiple stable 
modes. If the output converges to the target mode, 
the good response will keep the control parameter 
at /j. 0 -, if the output falls into some other mode, the 
change of the external response will automatically 
reset the control parameter to be /j c . Chaotic search 
can also be performed using the external reflectivity 
as the control parameter as shown in Fig. 6. Here, 
again three states (modes 9, 11, and 13) were set 
as targets. We purposely chose mode 11 as one of 
the targets to show a case when the search fails to 
converge. Since mode 11 is not among the mul¬ 
tistable states at fjP (see Fig. 2), the laser output 
cannot converge to this mode. This search method 
can only converge to modes which are stable at /x°. 

The influence of the parameter value for the 
onset of chaos on the search process has been in¬ 
vestigated. Figure 7 shows the effects of /x c on the 
average search time which was calculated from 50 
search events. It can be seen from Fig. 7 that the 
mode search time varies in a complicated way on 
the control parameter, the reflectivity. The mode 
search time depends on the details of the chaotic 
mode transition dynamics [Davis, 1998], which can 
have complicated structure depending on the con¬ 
trol parameter. More detailed analysis of this is a 
subject for future work. However, it is important 
to note that Fig. 7 shows that there is a parameter 
interval over which the variation of the search time 
is small enough for the search method to be con¬ 
sidered robust. On the other hand, for fi c < 2.1%, 
the average search time shows significant increase 
and the search fails when // < 2%. The increase 
of the search time near the boundary of the chaotic 
mode-transition regime, where the mode-transition 
rate decreases, shows that chaotic mode-transition 
is essential in the adaptive mode selection. 

We briefly discuss the feasibility of experi¬ 
mentally implementing the above mode selection 
method. For experimental convenience, we could 
substitute the mode-evaluation method used in the 



External Reflectivity /x c (%) 

Fig. 7. The dependence of average mode search time on 
control parameter p c . The control parameter is external 
reflectivity R ex t- Other parameters are h = 1.5 Ith and 
Text = 15 cm. Dependence on bias injection is similar. 

current simulation with a more simple method: To 
detect the light power for a particular wavelength 
and compare it with a preset threshold, i.e. to 
evaluate MPR instead of DMR. We numerically 
tested this simplified version of the search technique 
and obtained successful mode selection provided the 
threshold is appropriately chosen. The experimen¬ 
tal implementation could be realized by using a 
diffraction grating and a photodiode whose band¬ 
width corresponds to the mode-evaluation time, 
and an amplifier with a preset threshold voltage. 
The target mode can be set by varying the angle of 
the grating and hence the wavelength of the light 
incident upon the detector. The detected signal for 
a specific wavelength could be compared with the 
preset threshold and the deviation signal employed 
to adjust the injection current of the laser or the 
electric port of an acousto-optical modulator which 
controls the external reflectivity. 

4. Conclusion 

In summary, we have investigated mode-compe¬ 
tition dynamics of a multimode FP semiconduc¬ 
tor laser and demonstrated that different mode- 
competition states exist when external feedback or 
bias injections are varied and these states and their 
bifurcation structure are suitable for performing 
adaptive mode (wavelength) selection using chaotic 
search. For a solitary multimode laser, all modes 
are excited and no particular mode dominates the 
light output. When external feedback or injec¬ 
tion modulation are introduced, the laser exhibits 
complex mode-competition dynamics. Three differ¬ 
ent parameter regimes have been verified namely a 
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single-mode regime, a multistability regime, and a 
chaotic mode-transition regime. 

Switching between the multistability regime 
and the chaotic mode-transition regime makes it 
possible to adaptively select a dominant lasing 
mode by using a chaotic search method. In this 
method, the laser is initially set in the chaotic mode- 
transition regime. When a particular wavelength 
gets a good external response, the external response 
stabilizes that mode as the dominant mode of the 
laser by driving the control parameter (reflectivity 
or injection) into the multistability regime. We have 
performed computer simulations of the search tech¬ 
nique and demonstrated that the proposed method 
is quite effective over a wide parameter range in 
spite of its simplicity. It is noted that the proposed 
method may also be applicable to the strongly mul¬ 
timode class-B lasers [Szwaj et al., 1996] where 
hundreds of lasing modes oscillate simultaneously. 
Adaptive mode selection may provide a promising 
way for applying such lasers in adaptive communi¬ 
cation or adaptive data transmission networks. 
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We consider the effects of parametric perturbation on the onset of chaos in different dynamical 
systems. Favoring or suppression of chaos was observed depending on the phase or the frequency 
of the periodic perturbation. A lowering of the threshold of chaos was observed in an electronic 
device simulating a Josephson-Junction model and the suppression of chaos was obtained in a 
bistable mechanical device. We showed that in case of spatial instability in a sample of liquid 
crystal, the action of the parametric perturbation is to modify the velocity and the onset of the 
defects. 

Considering that the emergence of defects precedes the threshold of spatio-temporal chaos, 
we infer that parametric perturbation can modify the threshold of chaos in this spatial dynam¬ 
ical system. 


1. Introduction 

In 1965, an interesting paper appeared where 
Wehrmann [1965] was able to suppress turbulence 
behind a cylinder in a moving fluid. The basic idea 
was to put in vibration the cylinder with a suit¬ 
able feed-back using the same fluctuations present 
in the turbulent fluid. A complete laminarization 
was obtained. Turbulence is a phenomena related 
to a system with infinite degrees of freedom and 
it is natural to wonder if parametric perturbations 
can modify the onset of chaos in low dimensional 
system as well. In 1990, Lima and Pettini [1990] 
showed with rigorous theoretical consideration that 
resonant parametric perturbation can remove chaos 
in low dimensional systems. They confirmed this 
prediction with numerical simulation. Further¬ 
more, Cicogna [1990] showed, using a Melnikov 
Integral [Melnikov, 1963], how to modify the thresh¬ 
old of chaos by resonant parametric modulation. 


Considering that these works were essentially 
of theoretical nature, it was particularly in¬ 
teresting to try to verify these predictions in 
real systems with an experimental point of 
view. In spite of the development of other 
methodologies for the control of chaos, e.g. [Ott 
et al., 1990; Pyragas, 1992] and more recently 
[Boccaletti & Arecchi, 1995], the techniques for 
controlling chaos with parametric perturbations re¬ 
main a good method when it becomes impossible 
to apply feedback on the systems. In the next 
sections we will reassume the experimental results 
that have allowed us to verify these predictions 
[Cicogna & Fronzoni, 1990; Fronzoni et al, 1991]. 
The validity of this kind of methodology has been 
verified in other experiments performed in several 
laboratories [Azevedo & Rezende, 1991; Braiman k 
Goldhirsch, 1991; Meucci et al., 1994; Ding et al, 
1994; Chizhevsky k Glorieux, 1995; Chizhevsky k 
Corbalan, 1996]. 
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2. The Influence of Parameter 
Perturbations in a 
Josephson-Junction Model 

We considered a model that describes a wide class 
of phenomena that can be reduced to a pendulum in 
the presence of forcing and viscosity. This equation 
can be resumed in the following equation: 

+ S~- + D sin (f> + 7 cosiuit) = 0 (1) 

dt z at 

where 5 means the viscosity coefficient, 7 the am¬ 
plitude of forcing and uj the frequency. 

This system behaves chaotic in a range of 
parameters 7 and D as widely described in an 
experiment [D’Humieres et, al., 1982], where an elec¬ 
tronic device was used to simulate Eq. (1). 

We repeated this experiment but with an 
important improvement, which is to consider the 
coefficient D as a modulated control parameter: 

D = D 0 [l + £ cos(fit + 0)] ( 2 ) 



Fig. 1. Scheme of circuit for simulating Eq. (1). 


threshold of chaos. In this experiment we used forc¬ 
ing and modulation locked with 9 = 0. 

The prediction of the chaos threshold can be 
obtained by means of the Melnikov Integral. After 
some approximation [Cicogna & Fronzoni, 1990] one 
deduces a simple relation for this quantity 


where II means the frequency modulation, £ is the 
amplitude and 6 is the difference of phase between 
the modulation and the forcing. 

In Fig. 1 we show the scheme of the circuit 
which simulated Eq. (1) using minimum compo¬ 
nents technique [Fronzoni, 1989]. 

The main element of this device is the Voltage 
Control Oscillator (VCO) that provides an output 
with frequency depending on the input Vi according 
to the relation 


M(to) = —A cos(Qto) + B sen(Qfo + 6 ), (5) 

with 

A = 85 ± 27T7 sech(unr/2) ( 6 ) 

and 

B = 27r£fi 2 csch(fbr/2). (7) 

Assuming £ > 0, chaos does not appear if 

7 sech(cu 7 r/ 2 ) 4- £fi 2 csch(fl 7 r/ 2 ) < 4<5/ 7 T. ( 8 ) 


a; vco = KV i • (3) 

Ml and M2 are multipliers, and without going into 
great detail, it is important to know that the phase 
<p of the voltage at the VCO-output is described by 
the following relation: 


Vo(l + £ cos flf) S p - + 


(f> Vy C •• „ 

KRi + R^ + K (>) ~ ' 


(4) 


With suitable time-scaling this relation is equiva¬ 
lent to Eq. (1) which includes the modulation term. 
To study the influence of the modulation we fixed 
the amplitude and the frequency of the modulation, 
then the forcing was increased until it reached the 


It is important to observe that the sign of M(to) 
depends on the sign of £ and 9. In other words the 
threshold of chaos is defined by the relation of sign 
between the amplitude modulation £ and the phase 
9. For instance, £ > 0 and 9 > 0 provide a lowering 
of the threshold of chaos. 

Figure 2 shows the results, which are justified 
by this theory where one has to expect a lowering of 
the threshold of chaos versus the amplitude of the 
modulation. The quantitative differences between 
the experimental data and the theoretical predic¬ 
tions are due to the approximations used in the 
computation of the Melnikov Integral. 

When we performed this experiment we were 
not interested in observing the suppression of chaos, 
nevertheless the Cicogna’s theory predicts this pos¬ 
sibility. Only favoring chaos was experimentally 
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Fig. 2. Experimental results and theoretical predictions 
(z-line) for threshold chaos as function of the amplitude mod¬ 
ulation with 6 = 0, S = 0.25, to = 0.75 and fi = 0.75. 

verified with this kind of apparatus. Prodded by the 
theory of Lima and Pettini we performed another 
experiment with the specific intention to verify the 
suppression of chaos in a real mechanical system. 


3. Suppression of Chaos in a 
Bistable Mechanical Device 

The motivation behind this choice arises from the 
fact that the theoretical predictions could be in¬ 
validated by the incomplete correspondence of the 
model with the reality. This could cause doubts 
on the true possibility to remove chaos in a real 



Fig. 3. Scheme of mechanical device and the experimental 
configuration. 


system. For this purpose we restricted our choice to 
a bistable system as illustrated in Fig. 3 and mod¬ 
eled by a Duffing-Holmes forced oscillator [Moon, 
1997]. 

A metal beam was fixed over an oscillating 
plate. Two electromagnets are located near the 
free-end of the beam. This system is approximately 
described by this equation: 

4> = acj) — (3$ — 7 4> + A cos(27 ruft) . (9) 

The parameters a and (3 are obtained by a di¬ 
rect measurement of the equilibrium position ±xq 
and resonant frequencies uq for small perturbations, 
according to the following relations: 


a= i(27n/ 0 ) 2 , 

(10) 

/3= 2^ (27ri/ ° )2 - 

(11) 

We get the modulation by sending a periodic 
current into the electromagnets and this modifies 
the parameters a and f3: 

a —>• a[l + e cos(27a/ m t)], 

(12) 

/? —> f3[ 1 + 77 cos(27 TV m t)\ . 

(13) 


e and 77 are the modulation amplitudes and v m the 
frequency modulation. Figure 3 shows the experi¬ 
mental configuration used to control chaos. 

The shaker drove the plate with a frequency 
forcing Vf. An optical detector read the oscillation 
amplitude of the beam. This signal and its deriva¬ 
tive were sent to the x-y input of an oscilloscope. 
If the z-axis of the oscilloscope was triggered by the 
forcing, a Poincare section appeared on the screen. 
With this technique we were able to know in real 
time the state of the system without computations. 
Fixed dots on the screen indicated order and spread 
points on the screen indicated chaotic dynamics. 
According to the Lima-Pettini theory it is possi¬ 
ble to suppress chaos if the frequency of the pertur¬ 
bation approaches the values of the characteristic 
frequencies of the system, including harmonics. 

To observe the suppression of chaos we used 
the following procedure: We fixed the driving am¬ 
plitude and its frequency to observe a strange at¬ 
tractor on the screen and then we sent a periodic 
current into the coils. The frequency of this pertur¬ 
bation was swept around the values of the driving 
frequency or around the harmonics. 
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SECONDS 


Fig. 4. Time evolution of voltage at the output of Optical 
Detector in presence of modulation. Vf — 13.03 Hz, e = 
0.1, r] = 0.08, and A/A^ &os = 1.13. (a) v m = 13.08 Hz, 
(b) i'm = 13.07 Hz, (c) v m = 13.06 Hz. 


Figure 4 shows the voltage of the optical detec¬ 
tor versus the time for different modulation frequen¬ 
cies. The beam angle amplitude <p is proportional 
to this quantity. It results in: 


(a) Chaotic vibrations of the beam were alternated 
with ordered oscillations. 

(b) The ordered states were characterized by 
periodic oscillation around one of the two equi¬ 
librium positions. The beam amplitude was 
synchronized with the forcing. 


14 15 



28 30 V M (Hz) 


Fig. 5. Amplitude Lock-in output versus frequency mod¬ 
ulation I'm (a) first-harmonic resonance, upper curve e = 
0.16, lower curve e — 0.1. (b) second-harmonic resonance 
£ = 0.12. Sweep time 1175 sec/Hz, integration time 300 sec 
and AJAc haos 1.33. 


(c) The average duration of the ordered phases was 
proportional to the reciprocal of the difference 
Al/ = !//- u m . 

A way to get numerical evaluation of the res¬ 
onance phenomena, was to send the derivative of 
(^-voltage to the input of a Lock-in instrument us¬ 
ing the amplitude modulation as reference signal. 
Peaks in the response of the Lock-in indicate or¬ 
dered state synchronized with the modulation sig¬ 
nal. Figure 5 shows the Lock-in output versus the 
frequencies modulation. Resonance at the forcing 
frequency and at the second harmonics are evident. 
These results are well explained by the theory [Lima 
& Pettini, 1990]. 

4. Parametric Perturbation on 
Spatio-Temporal Instability 

The onset of chaos in spatio-temporal dynamics 
is characterized by the appearance of defects or 
topological defects. The behavior of these defects 
become more complex when a control parameter 
approaches the chaos threshold. For instance, the 
number and the velocity of these objects are increas¬ 
ing with the control parameter. 
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Due to the absence of regularity in the behav¬ 
ior of these defects its results are too ambiguous to 
define the threshold of the spatio-temporal chaos. 
In fact, in spite of the presence of a few defects in 
the structure, they assume chaotic motion just for 
a small increase of the control parameter. 

To study the influence of parametric perturba¬ 
tions on this kind of system, we considered an elec¬ 
trohydrodynamic instability [Williams, 1963] that 
appears in samples of nematic liquid crystals in 
the presence of electric field. Nematics are organic 
elements characterized by molecules with rod like 
shape. These properties induce an ordered configu¬ 
ration where the molecules prefer to stay in a well 
defined orientation around a mean value. A vec¬ 
tor n called “director” defines this orientation. In 
our experiment the nematic was closed between two 
conducting glass plates. The internal plate surface 
was rubbed in a well defined direction. In this con¬ 
figuration the molecules assume an orientation par¬ 
allel to the rubbing direction. When an electric field 
is applied on the sample, fluctuations of the molec¬ 
ular orientation induce spatial charges inside the 
nematic. The charges are generated by a current 
perpendicular to the field because of the anisotropy 
of the nematic conductibility. The interaction of 
the spatial charges with an electric field cause move¬ 
ment in the fluid. This motion arises over a critical 
threshold of the field and a regular pattern appears 



Fig. 6. Microscopic image of liquid crystals cell in pres¬ 
ence of Williams’ domains. Two defects are present in the 
structure. 



Fig. 7. Velocity of defects versus frequency modulation. 
Open circles correspond to the data of repeated experiment 
30 days later. 


on the sample, as shown in Fig. 6. The image is 
obtained observing the cell with a polarized micro¬ 
scope. Dark and white lines correspond to opposite 
velocities of the fluid in the sample with rotation 
axes parallel to the plates. As the electric field is 
increased topological defects broke the regularity of 
the structure. Two examples of these defects are 
shown in Fig. 6. These defects are well described by 
the Landau-Ginsburg theory [Ginsburg & Pitaevs, 
1958], 

In our experiment we used MBBA ( Methoxy- 
benzilidenebutylaniline) inside square glass plates 
of 2 cm size and separated by 20 ptm. We ap¬ 
plied on the sample an alternating voltage at the 
frequency of 10 Hz. We perturbed the system by 
means of voltage amplitude modulation of the order 
of 10%. Figure 7 shows the velocity of the defects 
against the frequency of periodic amplitude mod¬ 
ulation and a resonance behavior is evident. This 
experiment was repeated with the same sample 30 
days later obtaining a good reproducibility of the 
phenomena. 

It is important to note that a modulation of 
10% of the forcing at u m = 0.5 Hz induced an 
increase of the order of 50% of the defects veloc¬ 
ity. Considering that the increase of defects velocity 
precedes the onset of chaos means that parametric 
perturbation induces a lowering of the threshold of 
chaos in this kind of instability. 
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5. Conclusions 

In these experiments we obtained suppression or 
favoring of chaos depending on the frequency and 
phase of periodic perturbation applied on differ¬ 
ent systems. Moreover an interesting result was 
obtained in a sample of liquid crystal where the 
influence of perturbation modifies the velocity of 
defects in spatial structures. This behavior suggests 
that parametric modulation can modify the onset 
of chaos as a consequence of modifications of the 
dynamic properties of the defects in proximity to 
the turbulence threshold. In order to explain these 
last results we have developed a phenomenological 
model and are in the progress of further laboratory 
experimentation. 
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The paper considers the problem of designing time delayed feedback controllers to stabilize 
unstable periodic orbits of a class of sinusoidally forced nonlinear systems. This problem is 
formulated as an absolute stability problem of a linear periodic feedback system, in order to 
employ the well-known circle criterion. In particular, once a single test is verified by an unstable 
periodic orbit of the chaotic system, our approach directly provides a procedure for designing 
the optimal stabilizing controller, i.e. the one ensuring the largest obtainable stability bounds. 
Even if the circle criterion provides a sufficient condition for stability and therefore the obtained 
stability bounds are conservative in nature, several examples, as the one presented in this paper, 
show that the performance of the designed controller is quite satisfactory in comparison with 
different approaches. 


1. Introduction 

The problem of controlling chaos has recently re¬ 
ceived a lot of attention and many theoretical as 
well as experimental contributions coming from sev¬ 
eral areas have considered different aspects of this 
problem (see e.g. [Shinbrot et al., 1993; Ogorzalek, 
1993; Chen & Dong, 1993; Abed &: Wang, 1995; 
Hu et al., 1995] and references therein). One of 
the most frequent objectives consists in the sta¬ 
bilization of chaotic behaviors to periodic regimes 
and, in particular, many methods are directed to¬ 
wards the stabilization of one of the infinite un¬ 
stable periodic orbits that coexist in the chaotic 
attractor. This problem was first considered by 
Ott, Grebogi and Yorke in a paper [Ott et al., 
1990] that has originated, with related extensions 
and modifications (see e.g. [Romeiras et al., 1992; 
Hunt, 1991; Dressier &; Nitsche, 1992]), the so- 
called OGY methods. A distinct technique was pro¬ 
posed by Pyragas [1992], again forming the basis 


for other related works [Socolar et al., 1994; 
Kittel et al., 1995; Just et al., 1997], which are re¬ 
ferred to as time-delayed auto-synchronization or 
time-delayed or simply Pyragas methods. These 
two classical approaches, based on quite different 
properties and design procedures and successfully 
applied to various physical experiments, lead to 
feedback controllers utilizing small perturbation 
signals and they really exploit the peculiar char¬ 
acteristics of chaos. Also for this reason, the above 
methods and in general most literature on control of 
chaos appear to be distinct from what is called con¬ 
trol theory, although some evident links are present. 

The purpose of the present paper is to 
contribute (see also [Basso et al., 1997a]) to the 
results of feedback control theory in order to opti¬ 
mize, in a certain sense, the solution obtainable via 
the Pyragas approach, and also to indicate a means 
of designing the related controller. In particular, 
we consider periodically forced chaotic systems, for 
which the Pyragas methods are quite appealing, 
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and we formulate the periodic orbit stabilization 
as an absolute stability problem [Vidyasagar, 1992] 
of a linear periodic feedback system. Therefore, 
a well-known frequency criterion, the circle crite¬ 
rion [Vidyasagar, 1992], is employed to compute 
the largest obtainable stability bounds. Once a sin¬ 
gle test is verified by the unstable periodic orbit 
of the chaotic system, our approach directly pro¬ 
vides a procedure for designing the “optimal” sta¬ 
bilizing controller, i.e. the one maximizing the sta¬ 
bility bounds. Even if the method is based on a 
sufficient condition and therefore the bounds are 
conservative in nature, several examples show that 
the performance of the optimal controller is quite 
satisfactory. 

The remainder of the paper is organized as fol¬ 
lows. Section 2 formulates the problem and de¬ 
velops the main results, while Sec. 3 presents in 
some detail one example showing the application 
of these results, also in comparison with different 
approaches. Some brief comments end the paper in 
Sec. 4. 


2. Main Results 

The notation adopted hereafter is pretty standard. 
For instance, X(s) denotes the Laplace transform 
of the time-domain signal x(t), f stands for the 
derivative of a function /. A transfer function is 
said to be stable if its impulse response belongs to 
Li, the space of time-domain signal x(t ) such that 
/ 0 °° \x(t)\dt is finite. If the transfer function is ra¬ 
tional, this reduces to the usual notion of stability, 
i.e. all the poles must have a negative real part. 

To motivate our problem of stabilizing unsta¬ 
ble periodic solutions, we consider one of the most 
studied problems in controlling chaos, i.e. the elim¬ 
ination of a period doubling bifurcation sequence 
[Abed & Wang, 1995]. 

Consider the periodically forced (Lur’e) sys¬ 
tem depicted in Fig. 1, where G(s ) is the transfer 
function of a stable finite dimensional linear time- 
invariant system and n(-) is a sufficiently smooth 
nonlinear function such that n(0) = 0 and n'(0) = 
0. These latter hypotheses on the function n(-) 
are quite general (see also the example of Sec. 3) 
and are posed to ensure that the origin is a locally 
stable equilibrium point of the uncontrolled system 
(i.e. u(t) = 0), when the amplitude A of the forcing 
term is zero. For increasing values of the ampli¬ 
tude A the uncontrolled system exhibits a stable 


A sin 


+ 

s 

G{s) 


) - 






n(-) 





Fig. 1. Feedback system. 



Fig. 2. The solid (dotted) curves correspond to stable (un¬ 
stable) T-periodic solutions. (t/max := max (s [o,r)3M(t), 
1/min := min( S [ 0 j T] IM(*))- 


T-periodic solution up to a certain value, where a 
first period doubling bifurcation appears. After this 
value the system possesses an unstable T-periodic 
solution and a stable 2T-periodic solution, until this 
latter solution undergoes a new period doubling bi¬ 
furcation, and so on according to the well-known se¬ 
quence which usually leads to chaos. This scenario 
is illustrated in Fig. 2 where only the T-periodic so¬ 
lutions are indicated explicitly. The basic idea for 
eliminating the bifurcation sequence is to stabilize 
the unstable T-periodic solution (the dashed one 
in Fig. 2), thus inhibiting the first period doubling 
bifurcation. Obviously, this calls for the design of 
a controller that stabilizes the family of T-periodic 
solutions corresponding to the largest possible in¬ 
terval of amplitudes A of the forcing term. 

To proceed, consider the problem of stabiliz¬ 
ing an unstable T-periodic solution j/a( 0 of the 
uncontrolled system (i.e. u(t) = 0) of Fig. 1 cor¬ 
responding to a given amplitude A of the forcing 
term. Such a problem requires the design of a linear 
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time-invariant feedback controller of transfer func¬ 
tion C(s ), i.e. 

U(s) = C(s)Y(s), (1) 

such that yA(t) becomes a stable T-periodic solu¬ 
tion of the controlled system. 

It is easily recognized that the controller trans¬ 
fer function must satisfy the condition 

C (jn-fr) =0 f° r n = °> 2 ’'' • ( 2 ) 

to ensure that the T-periodic solution y A (t ) is still a 
solution of the controlled system. This fact implies 
that the transfer function of such a controller has 
the form 


It is well known that the stability of y A (t) is guar¬ 
anteed by the stability of the origin of the linearized 
system of Fig. 3 [Vidyasagar, 1992]. Here, we refer 
essentially to input-output stability, even if quite 
mild assumptions on the feedback interconnection 
of G(s ) and C(s) allow one to consider asymptotic 
stability [Basso et al., 1997a]. 

To assess the stability of the system of Fig. 3 
and consequently the stability of y^(£), we can em¬ 
ploy the well-known circle criterion [Vidyasagar, 
1992]. The two quantities 

1 1 

a :=--——r p :=-;--—— 

maxjqo , T }k A {t) mm te [ 0 j t] (t ) 

( 6 ) 

are introduced to state the following stability result. 


C(s) = F(«)(l - e- T ) (3) 

where F(s) is in general not restricted to be ra¬ 
tional. For instance, F(s) is a constant gain K 
in [Pyragas, 1992], while it has the form K{ 1 — 
Re~ sT )~ l in [Socolar et al., 1994], where R is a 
suitable constant. A generalization of this latter 
structure has been used in [Basso et al., 1997a]. 

To investigate the stability of y A {t), the lin¬ 
earization of the Lur’e system around the periodic 
solution is performed. This leads to the linear pe¬ 
riodic system of Fig. 3, where the transfer function 
L(s ) of the linear subsystem is given by the feedback 
interconnection of G(s) and C(s), i.e. 


Theorem 1. Let a < 0 and (3 > 0. Then, the con¬ 
troller C(s) stabilizes yA(t) if 

1. L(s) is stable; 

2. the following inequality holds: 


L(ju>) - 


a + (3 
2 


< 


(3 — a 
2 


Mu > 0. (7) 


Remark 1. Condition (7) has a simple graphical in¬ 
terpretation in terms of the Nyquist plot of L(s). 
Indeed, the Nyquist plot must lie inside the circle 
of center ((/3 + a)/ 2, 0) and radius {(3 — a)/ 2. 


L(s) 


G(s) 

1 + C(s)G(s) 


(4) 


and the periodic gain in the feedback path has the 
form 


Remark 2. Similar results can also be given for gen¬ 
eral values of a and f3. However, for the purposes of 
this paper, it is enough to consider the case a < 0 
and (3 > 0. 


, . . dn 

hA{t) = Ty 


y=VA(t) 


(5) 



From condition (2) and Eq. (4), it follows that 
L(jn2n/T ) = G{jn2n/T ) for n = 0,1,2,.... 
Therefore, it is evident that Theorem 1 can ensure 
stability of y A {t) only if L(s) is stable and the fol¬ 
lowing inequalities 


G 



a + (3 
2 


B — OL 

< —-— for n = 0, 1,2,- 

Z 

( 8 ) 


hold. Note that the inequalities above rely on the 
uncontrolled system, thus resulting in an a priori 
test for the guaranteed stabilizability of y^i(t). 

Now, the major question is the following: If the 
a priori test (8) is satisfied for yyi(t), is it possible 
to design a controller C{s ) such that Theorem 1 
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ensures stability of t/^(f)? Fortunately, the answer 
is affirmative as shown by the next result. 


Theorem 2. Let pa(1 ) satisfy condition (8) and 
define $(s) = 1 — e~ sT . Then, it is always pos¬ 
sible to find a stable rational (with stable inverse ) 
approximation $ ap (s) of$(s) and a stable rational 
transfer function W ( s ) such that the controller 


C(s) = 


1 - 


a + /3\ $(s) 


G(s) 


2 G(s)J $ ap (s) 


IT(s) 


1 - 1 - 


a + f3\ $(s) 


(9) 


2 G(s)J $ ap (s) 


W(s) 


stabilizes y^(t). 

Sketch of the proof (see [Basso et al., 1997b] for 
details): 

The stability of L(s) follows directly from the 
stability of G(s), $~ p (s), W(s) and 3>(s). Indeed, 
from Eqs. (4) and (9) it turns out that 

L W = G(S )-( GW -^)^ W . 

Consider the inequality (7). After some manipula¬ 
tions, it can be rewritten in the form 




\ju) 
Vw > 0 


f3 - a 


and, equivalently, 




W (ju) 


j3 — a 


< 


G(ju) - 


a + (3 


Mu > 0. 


( 10 ) 


Now, the ratio at the right hand side of inequality 
(10) is greater than one for the set of u such that 
G(ju ) is inside the circle of center ((/?+a)/2, 0) and 
radius ({3 —a )/2 (see Remark 1). For this set, which 
contains u — n2n/T for n = 0, 1, 2,..., as ensured 
by condition (8), it is enough to put W (ju) close 
to zero in order that condition (10) holds. For the 
remaining values of u, i.e. those such that the ra¬ 
tio at the right hand side of inequality (10) is equal 
or less than one, we select 4> ap such that the ra¬ 
tio $(ju)/<& ap (ju) is arbitrarily close to one. This 
can always be done since at such frequencies <&(ju) 


is always different from zero. At this point it is 
enough to choose W (ju) close to one to fulfill con¬ 
dition (10) also for these frequencies. The fact that 
rational transfer functions are dense in the L\ space 
concludes the proof. 

Remark 3. A possible way to select 4> ap is to use the 
standard odd Pade approximation of e~ sT . How¬ 
ever, the resulting $ ap has its zeros on the imagi¬ 
nary axis. This drawback is simply eliminated by 
slightly shifting the zeros into the left half plane, 
thus arriving at a new $ ap that has a stable inverse. 
Increasing the order of the Pade approximation, it is 
easily verified that the ratio <&(ju)/<b- dp (ju) is quite 
close to one at all frequencies, except for u = n27r/T 
for n = 0, 1, 2,.... The choice for W is based on the 
idea of not using the controller (W(ju) — 0) at the 
frequencies u where G(ju) is already inside the cir¬ 
cle of center ((/3+a)/ 2, 0) and radius (fi—a)/ 2, and 
to use the controller at the remaining frequencies 
to force L(ju) towards the center of such a circle. 
Finally, we note that it is often possible to simplify 
the controller (9) via a suitable order reduction. 

Despite the complicated aspect of Eq. (9), the 
implementation of the controller is rather simple. 
Indeed, the controller is obtained via the positive 
feedback scheme reported in Fig. 4 that makes clear 
the presence of a unique delay element and two ra¬ 
tional filters, G(s) and 

/ a + <n rn«) 

’ \ 201.,)) G(s)<b ,„(>) ■ 

Let us now consider the problem of stabilizing 
the family of T-periodic solutions corresponding to 
the largest interval of amplitudes A of the forcing 
term. To this purpose, let Aq denote the largest 
possible value of A such that condition (8) holds 
where a = ao and — /?o, being ao and fio defined 
as 

ao = max a fio = min (3. 

A^[0,Ao] Ae[0,.Ao] 





Y{>) 


Fig. 4. Controller structure. 
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Clearly, the application of Theorem 2 guarantees 
that all the T-periodic solutions corresponding to 
amplitudes A < Aq are stabilized by the controller 
C(s). 

Finally, we recall that Theorem 2 provides 
a sufficient condition for the stability of yA(t), 
i.e. i/a(£) may be stable even if condition (8) is far 
from being satisfied (see the example in the next 
section). However, several application examples 
have shown that the controller designed via our ap¬ 
proach performs much better than other controllers 
that do not satisfy condition (7). This is illustrated 
via the example developed in the next section. 


3. Example 

Consider the equation of the driven Toda 
oscillator 


y + 0.8 y + 25(e v — 1) = A sin 2nt — u(t). (13) 


This system has been extensively studied in [Just 
et ai, 1997], where the analysis of the Floquet ex¬ 
ponents showed the presence of a set of stable orbits 
of period T = 1 in the range A € (0, 37) of the forc¬ 
ing input. For larger values of the amplitude the 
system undergoes a sequence of period doubling bi¬ 
furcations leading to chaos. The stable (solid curve) 
and unstable (dotted curve) T-periodic solutions of 
the uncontrolled system are reported in Fig. 2. The 
upper and the lower curves represent the maximum 
and minimum, respectively, of the periodic solution 
2 M(t), as a function of the system parameter A. 

The control technique proposed in this paper 
can be fruitfully exploited to derive a unique time- 
delayed feedback controller such that the controlled 
system has stable periodic solutions for a much 
larger range of the parameter A. It is easily recog¬ 
nized that system (13) can be recast in the classical 
Lur’e structure of Fig. 1, where 


G(s) 


1 

s 2 + 0.8s + 25 


(14) 


and 


n(y) = 25(e y - y - 1). (15) 



Fig. 5. ot(A) (lower curve), /3(A) (upper curve). 


0.25 
0.2 
0.15 
0.1 
0.05 
' 0 
- 0 . 051 


- 0.1 


- 0.15 



- 0 . 2 - 
- 0.25 - 

_i_i_i_i_i_i_i_ 

- 0.3 - 0.2 - 0.1 0 0.1 0.2 0.3 

Real 


Fig. 6. Nyquist plot of L(s) using controller (16). 


From the bifurcation diagram of Fig. 2 and us¬ 
ing expressions (5), (6) and (15) we can derive the 
two functions a and (3 of the amplitude A, shown 
in Fig. 5. 

From these functions, we compute Ao, i.e. the 
largest possible value of A satisfying condition (8), 
that is the smallest circle that includes the fixed 
complex points at the frequencies 2mr, n = 0, 1,... 
(‘x’ marks in Figs. 6-8). It turns out that «o = 
-0.064, /3b = 0.088 and A 0 = 8.6. 

Based on Theorem 2 and the structure in Fig. 4, 
the following controller is obtained 


C(s) 


H(s)( 1-e-fy 
1 - G(s)H(s){l - e~ s ) ’ 


Note that G(s) is stable and n(0) = n'(0) = 0. 


(16) 
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where 


H(s) = 


1.3s 2 + s + 32.5 


5 • 10' 7 s 3 + 1.3 • 10~ 4 s 2 + 1.6 • 10 _2 s + 1 

■*£(*)■ (17) 

Here, the stable filter 


^ap( 5 ) 

2s 5 + 20s 4 + 899s 3 + 4598s 2 + 35720s + 108900 
~ s 5 + 30s 4 + 420s 3 + 3360s 2 + 15120s + 30240 

(18) 

follows directly from the fifth order Pade approxi¬ 
mation of e~ s . 

The Nyquist plot of the corresponding closed 
loop transfer function L(s) is indicated in Fig. 6 to¬ 
gether with the circle intersecting the real axis in 
ao and /? 0 . 

Simulations show that the designed controller 
(16) stabilizes the periodic orbits of system (13) 
over the range A € (0, 1500), where the upper 
bound is only due to numerical problems in inte¬ 
grating the system at large values of A. 

We can compare the results obtained by our 
controller with those presented in [Just et al. , 1997], 
where a Pyragas controller was designed exploiting 
Floquet theory. The derived control law was 

u{t) = 2A[y(t) - y(t - 1)], (19) 

or, equivalently, 

C{s) = F(s)( 1 - e~ s ) = 2.4s(l - e~ s ) . (20) 

Such a controller is optimal among the con¬ 
trollers with the same structure but different gains 
(i.e. F(s) = Ks ), in the sense of minimizing the 
Floquet exponents at the given input amplitude 
A = 105. 

Figure 7 shows the application of Theorem 1 to 
the system (13) and the controller (20). From the 
Nyquist diagram we can derive A = 4, a = —0.149 
and (3 = 0.173 (circle with solid curve). For com¬ 
parison, Fig. 7 also reports the optimal circle at 
Ao = 8.6 (dotted curve). 

Finally, we note that the controller (20) is also 
optimal with respect to Theorem 1 in the sense that 
maximizes the amplitude A among the controllers 
with the same structure but different gains. 

Simulations on the controlled system show that 
the controller (20) can stabilize all the orbits of pe¬ 
riod 1 in the much restricted range A 6 (0, 131). 



Fig. 7. Nyquist plot of L(s) using controller (20). 



Fig. 8. Nyquist plot of G(s). 


Table 1. Results of the example. 

Guaranteed Simulated 

Stability Stability 

C{s) Range Range 

0 A £ (0, 2.45) A £ (0, 37) 

2.4s(l - e~ s ) A £ (0, 4) A £ (0, 131) 
Eq. (16) A £ (0, 8.6) A £ (0, 1500) 


To conclude, Theorem 1 is applied to the un¬ 
controlled system [see the Nyquist plot of the linear 
subsystem (14) in Fig. 8] determining the values of 
a = —0.241, (3 = 0.266 and A = 2.45 (solid circle in 
Fig. 8). 

We summarize the results obtained in 
Table 1, which definitely supports our claim that 
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-40 -20 0 -40 -20 0 -40 -20 0 

(A=140) (A=140) (A=140) 



-400 -200 0 -400 -200 0 -400 -200 0 

(A=1200) (A=1200) (A=1200) 


Fig. 9. Simulations in the ( y , y) plane: Each row denotes a different value of the parameter A (40,140,1200), while each 
column indicates a different controller, i.e. column 1 C(s) — 0, column 2 controller (20), column 3 controller (16). 


the larger the amplitude of the system input fulfill¬ 
ing Theorem 1, the larger the range that is stabi¬ 
lized in practice. For completeness, we have also 
reported some of the simulation results in Fig. 9 
for different values of the parameter A and for the 
controllers presented in the example. 

4. Conclusion 

The paper has dealt with the problem of designing 
time delayed feedback controllers to stabilize un¬ 
stable periodic orbits, a problem of great interest 
in the control of chaos. For an important class of 
sinusoidally forced nonlinear systems, we have for¬ 
mulated such a problem as an absolute stability 
problem of a linear periodic feedback system, in or¬ 
der to employ classical stability tools. 

Our main contribution is a simple procedure 
for designing the optimal stabilizing linear time- 
invariant feedback controller, i.e. the one that guar¬ 
antees the largest stability bounds obtainable via 
the well-known circle criterion. Several examples 
have shown that the performance of the designed 
optimal controller is quite satisfactory also in com¬ 
parison with different approaches. This has been 
illustrated in the application example, where also 
the possible use of the proposed approach for elimi¬ 


nating the period doubling bifurcation sequence has 
been shown. 
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It has been established that chaotic iterated maps can be transferred to a periodic motion 
by application of periodic signals. As a next step, we investigate the optimization of such 
signals with respect to a small norm. This method employs a cost function on the space of 
periodic inputs. By means of a direct optimization search in this space, the approach stabilizes 
periodicity with small forces. The results for increasing driving periods are compared. 


1. Introduction 

In recent years, many authors considered the prob¬ 
lem of controlling chaotic systems (see [Chen &: 
Dong, 1993; Shinbrot, 1995] for reviews). A 
common goal is to force the chaotic system into a 
periodic state. Other intentions might be a migra¬ 
tion of the system to certain phase space regions, 
synchronization with an arbitrary signal, or the en¬ 
hancement of chaos in the system. In this paper 
we will focus on the conversion of chaos to peri¬ 
odicity. There are two aspects of this approach: 
(i) some advantage of a periodic behavior over a 
chaotic one (e.g. reduced drag, higher mean out¬ 
put, limited peak values, or simple predictability), 
and (ii) the proximity of the chaotic motion to 
(infinitely many) unstable periodic solutions. The 
second point is the reason why it may take very 
little effort to control the system. This paper inves¬ 
tigates a method to systematically reduce periodic 
control inputs that remove chaos in iterated maps. 

Removal of chaos with very little effort by use 
of its own specific complex phase space structure 
is the main aspect pursued in the chaos control 
field. In this regard, many proposed methods aim 
to exactly stabilize a pre-existing unstable periodic 
orbit (UPO) embedded in the chaotic attractor. 


Such techniques usually require some type of feed¬ 
back since the control forces correct deviations from 
the desired orbit. In principle (i.e. in absence of 
noise, modeling errors, etc.) control forces vanish 
if the goal is reached (“noninvasive” controls). By 
contrast there are “invasive” methods that always 
influence the system. The particular invasive tech¬ 
nique we consider is called “open-loop” since it fore¬ 
goes any additional feedback. The control forces 
are applied permanently and irrespective of the ac¬ 
tual system state. If one wants to achieve a peri¬ 
odic goal motion by an open-loop control, the forces 
have to be periodic, too. The only possible goal 
trajectories are asymptotically stable solutions of 
the periodically driven system — which are typi¬ 
cally not pre-existing unstable system solutions like 
UPOs. It might be possible, however, to control 
a goal close to a UPO. For instance, a pendulum 
can be stabilized in a vibrating state near to its un¬ 
stable inverted position by periodic up and down 
motion of the suspension (vibrational control, see 
e.g. [Bellman et al., 1986]). 

The history of periodic control of chaotic sys¬ 
tems dates back to Alekseev [1987] who sinusoidally 
modulated a parameter of a chaotic population 
system modeled by an ordinary differential equa¬ 
tion. More recent work has focused on driven 
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oscillators, often employing an analytical treat¬ 
ment by Melnikov’s method [Lima & Pettini, 1990; 
Braiman & Goldhirsch, 1991; Chacon, 1995]. Inter¬ 
esting effects of periodic and stochastic parameter 
modulation of iterated maps have been investigated 
by Rossler et al. [1989]. Subsequent experimental 
work has employed periodic modulation for removal 
of chaos [Azevedo & Rezende, 1990; Fronzoni et al ., 
1991; Ciofini et al ., 1995], 

A somewhat different approach of open-loop 
control was proposed by Hiibler k Liischer [1989]. 
They constructed control forces starting from an 
explicit goal trajectory of the system (whereas sim¬ 
ple periodic modulation methods do not provide a 
priori knowledge of the eventual controlled orbit). 
Their entrainment control method is not limited to 
chaotic systems, and the goal trajectory need not be 
periodic. Nevertheless, authors often deal with the 
chaos to periodicity conversion problem [Jackson k 
Hiibler, 1990; Jackson, 1991; Jackson Sz Kodogeor- 
giou, 1992; Mettin et al., 1995]. 

Modulation can be simple and fast, which is the 
main advantage of the open-loop approaches over 
the feedback approach as long as one accepts larger 
permanent control signals. An attempt to reduce 
the control power by a suitable shaping of the mod¬ 
ulation was proposed for chaotic continuous flows 
[Mettin k Kurz, 1995]. It was demonstrated that 
multimode signals can scale down the control and 
switch between different goal states. Now we con¬ 
sider the control of chaotic iterated maps in the 
same manner. Therefore, the control problem is 
restated in terms of an optimization problem. We 
seek to find convenient periodic control inputs that 
(i) introduce asymptotically stable periodic solu¬ 
tions in the chaotic system, and (ii) are small in 
a given sense. With respect to smallness, there 
are several considerations. In principle, the power 
needed by noninvasive methods to keep an origi¬ 
nally chaotic system in periodic motion is negligi¬ 
ble. Intuitively, one expects that the finite power 
an open-loop control needs is less the closer the 
goal is to the uncontrolled, natural dynamics. This 
might depend on the exact definitions of “power” 
and “close”. However, if this expectation holds, or¬ 
bits resulting from small signal periodic modulation 
may emerge near to UPOs. On the other hand, the 
underlying UPO structure may facilitate a periodic 
modulation approach. 

The optimization problem and the control sig¬ 
nal search are stated in detail in Sec. 2. In Sec. 3 
the method is applied to a map that models driven 


dissipative oscillators, and the results for increasing 
control signal period are compared. A discussion is 
given in Sec. 4. 

2. Optimization Problem and 

Direct Search Approach 

We consider an iterated map dynamics subject to a 
periodic control signal input, 

x (n) = f ( x (n) ) u (n)) 5 U (n+N) = u (n) ? 

where n <E Z is the discrete time, N is the period of 
the control, and x G u £ R fc , f : l d x I 1 ’ 4 R d . 
Without control, i.e. if all = 0, we assume 
chaotic motion of the system. The task is to find 

a convenient set U/V = { U W} n= o. N—l such that 

(i) there exists an asymptotically stable periodic so¬ 
lution of Eq. (1), and (ii) some norm of the input 
set is small. 

To determine the quality of a test set U/y, in¬ 
formation must first be gained about existence of 
an asymptotically stable periodic solution under its 
action. This is used in connection with the control 
signal’s norm for evaluation of a cost function like 

f ||Uiv|i : stable period 
cost = < (2) 

( c p : otherwise 

where Cp is some large (penalty) number. An op¬ 
timization algorithm performs the task of finding 
a test set that minimizes the cost function. It 
turns out that the cost landscape defined by Eq. (2) 
can be very complicated and rough or fractal-like. 
Therefore we expect to find sufficient but not opti¬ 
mal control forces. 

The optimization works directly on the ( kN )- 
dimensional space of control signal sequences Uat. 
The detection of the existence of an asymptoti¬ 
cally stable periodic solution is done by a recurrence 
check for different initial conditions. For simplicity 
in this initial study, we do not include statements 
about the basin of attraction so as to avoid excessive 
numerics needed to check for coexisting attractors. 
Thus we are content with periodicity for a few se¬ 
lected initial states xand an a posteriori test 
of the basin after the complete optimization. If all 
tested initial x( 0j - ) lead (after M transient itera¬ 
tions) to a periodic recurrence of the same response 
period N, it is concluded that a stable periodic solu¬ 
tion exists under the influence of the control signal, 
and that the basin of attraction has some larger ex¬ 
tent. Then the tested control signal is assigned a 
cost value according to its norm. 
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To decrease the cost, the optimization algo¬ 
rithm directly manipulates the control sequence U n 
what we call a direct search. It is also possible to 
implement an indirect search that manipulates goal 
trajectories instead of the control forces. The latter 
corresponds to a transfer of the optimized entrain¬ 
ment control approach [Mettin, 1997] to iterated 
maps, and it will be investigated in future work. 


3. Example 

As an example we consider a two-dimensional map 
introduced by Parlitz et al. [1991]. It models a class 
of driven dissipative oscillators and reads 

x[ n+1 ) = cos j3 — x fj 71 ) sin 0) 

^(n+i) _ s j n p _|_ x ( n ) cos ^3) + a + u 

where u^> is the (scalar) control, and 

7 = e~ dT ! 2 , /? = r+|4” ) l 2 + l4“ , l 2 - (4) 

The parameters correspond to a damping d, a driv¬ 
ing period T, and a driving amplitude a (where 
the control acts upon). They were set to T = 1.0, 
d = 0.9162, a = 1.5 which results in a chaotic dy¬ 
namics. As a cost function we chose Eq. (2) with 
Cp = 10 and the L 2 norm of the control inputs Uat 
written as an A'-dimensional vector, divided by the 
control period N: 

1 /' N - 1 \ 5 

ll u «H = ]v(X> < “ > l ! ) ■ (5) 

The recurrence was checked up to period 32 with 
five different initial conditions x(° J ) = (2 j — 1) 
(0.1,0.1), j = 1,..., 5. The check was done af¬ 
ter M = 100 transients and for a recurrence of 
10~ 3 . For the search of small cost values in the N- 
dimensional control signal space the numerical op¬ 
timization algorithm amebsa, a variant of simulated 
annealing, from [Press et al., 1992] was employed. 

Figure 1 shows a bifurcation diagram of the pa¬ 
rameter a near the chosen value for the uncontrolled 
oscillator map. It can be seen that the chaotic re¬ 
gion reaches from about a = 1.32 upwards, and 
a few, narrow periodic windows are visible. Peri¬ 
odic control is equivalent to a periodic shift of the 


i- 1 - 1 - 1 -r 



1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 


a 

Fig. 1. Bifurcation diagram of the uncontrolled oscillator 
map Eq. (3): first coordinate xi versus parameter a. 

driving amplitude a. Note, however, that a is never 
reduced below 1.32 by application of all in the fol¬ 
lowing reported optimized signals. 

A period 1 control signal corresponds to a con¬ 
stant parameter shift of a, and the nearest low- 
periodic window to a = 1.5 can be found at a\ = 
1.4356 (response period N = 28). This is in some 
sense the only obvious nonoptimized control signal 
as it can easily be extracted from a codimension 
one bifurcation diagram. Such an initial guess is 
not straightforward in the higher dimensional con¬ 
trol signal spaces (besides the trivial set of identical 
u = a - ai = 0.0644). Therefore the purpose of 
the optimzation procedure is twofold: first to find a 
nontrivial stabilizing sequence Ujv at all, and sec¬ 
ond to reduce the norm of such a sequence. 

Control inputs of periods N from 2 to 12 were 
used for optimization. For all control periods N the 
algorithm easily found nontrivial stabilizing sets. 
Further optimization reduced the norm on average 
by a factor of 2. The computations were terminated 
when no further improvement of the cost criterion 
(2) could be found within a given time. The runs 
were repeated a few times, and arbitrary sequences 
as well as results from lower periods (extended to 
an integer multiple control period) were used as ini¬ 
tial guesses. The overall best results from all runs 
have been taken. The expenditure for obtaining the 
individual optimized results was moderate but sim¬ 
ilar for the various N. Therefore a comparison of 
these is justified, though the results are supposed 
to provide only some local cost minima. Figure 2 
presents the magnitude of the best signals found in 
different norms versus the driving period N. The 
^ 2 -type norm used for optimization [Eq. (5)] is in¬ 
dicated by the squares connected by the solid curve. 
The dashed line connects for the same control in¬ 
puts the arithmetic mean of the absolute values 
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Fig. 2. Magnitude of the optimized control signals versus 
their period: the solid line connects the values of norm 
Eq. (5), the dashed line of norm Eq. (6). 


(Li-type norm): 

ll u »ll. = T? Z> < “ ) l- W 

71=0 

Note that Eq. (6) is a better criterion for compar¬ 
ison with the trivial control (N = 1) than Eq. (5) 
since the latter gives smaller values for the same 
signal written in a multiple period. The numbers 
above the squares indicate the actual system re¬ 
sponse period N that is achieved by control (always 
an integer multiple of the driving period N). All 
control signals have rather small norm values, and 
the resulting stable periodic orbits turn out to lie 
close to the chaotic attractor. It can be seen that 
there is the tendency of the norm to decrease for 
higher periods, which is plausible because of more 
degrees of freedom in the search space. However, 
the control cost does not decline monotonously, but 
seems to saturate except for some selected minima 
(at control periods N = 5, 7, and 10). Therefore the 
gain of an optimized periodic driving with respect 
to constant parameter shift (trivial periodic control) 
depends to some extent on the specific period of 
the applied control signal. In Fig. 3, the optimized 
period-7 control is shown (lower plot). The result¬ 
ing orbit was detected as a period 14, and its first 
coordinate is shown in the upper plot (diamonds 
connected by the solid line). It runs close to a 
UPO of period 7 (dotted line). In fact, all emerging 
stable periodic states are adjacent to pre-existing 
unstable periodic orbits, but a detailed analysis of 
this is beyond the scope of the paper. In particu¬ 
lar, a certain hierarchy of the different UPOs might 



Fig. 3. Optimized control signal u of period 7 versus discrete 
time n (lower plot), and resulting periodic orbit (diamonds in 
the upper plot). An adjacent UPO is indicated by the dotted 
line. 

underly the effort necessary to stabilize the system 
in their vicinity. 

Finally, the result of the optimization depends 
on the chosen cost function. Here a quadratic mean 
of the applied signal was minimized, which does not 
avoid relatively large single inputs (compare the rel¬ 
ative peak of u ^ shown in Fig. 3). To minimize 
the largest control amplitude, one might choose the 
maximum norm for the cost evaluation. 

4. Discussion 

It was demonstrated that a chaotic iterated map can 
be controlled to return eventual periodic motions by 
means of small optimized inputs of various period¬ 
icity. The optimization of the cost function corre¬ 
sponds to a search for periodic windows in the con¬ 
trol signal parameter space near to the origin. The 
possibility of finding such a window should in prin¬ 
ciple increase with the dimension (the control signal 
period). Accordingly, the example shows a trend of 
decreasing control forces for higher driving periods, 
with some periods being preferred. This might be 
related to the underlying structure of UPOs, as con¬ 
trolled trajectories appeared near them. Besides the 
pure reduction of open-loop control forces, this con¬ 
nection of stabilized periodic states and pre-existing 
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UPOs provides an interesting result. Further inves¬ 
tigations about the stabilizing mechanisms would 
be relevant not only for control applications but 
also for synchronization studies. Besides the role 
of unstable orbits, the investigation of chaotic tran¬ 
sients of the periodically controlled maps might be 
of some further interest since they can strongly vary 
for neighbored initial conditions. 
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We study the spatiotemporal dynamics of the underdamped Josephson junction series arrays 
(JJSA) which are globally coupled through a resistive shunting load and driven by an rf bias 
current. Clustering bifurcations are shown to appear. In particular, cluster-doubling induced 
period-doubling bifurcations and clustering induced spatiotemporal chaos are found. Further¬ 
more, an interesting spatiotemporal intermittency is also found. These phenomena are closely 
related to the dynamics of the single cell. 


The dynamics of globally chaotic systems has been 
of great interest in recent years. They arise natu¬ 
rally in studies of Josephson junctions arrays, mul¬ 
timode laser, charge-density wave, oscillatory neu¬ 
ronal system, and so on. Some rather surprising 
and novel features, such as clustering, splay state, 
collective behavior, and violation of the law of large 
numbers are revealed in these continuous and dis¬ 
crete globally coupled models [Benz et al. , 1990; 
Bhattacharya et al., 1987; Chernikov &; Schmidt, 
1995; Dominguez et al., 1991; Dominguez &; 
Cerdeira, 1995; Eikmans & van Himbergen, 1991; 
Fisher, 1983; Free et al., 1990; Hadley Sz Beasley, 
1987; Hadley et al., 1988; Hebboul &: Garland, 1991; 
Kaneko, 1989; Kvale & Hebboul, 1991; Lee et al., 
1992; Middleton et al., 1992; Strogatz & Mirollo, 
1993; Tchiastiakov, 1996; Tsang et al., 1991; Tsang 
&: Schwartz, 1992; Watanabe &; Strogatz, 1993; 
Wiesenfeld et al., 1996]. 

Being a paradigm for the study of nonlinear 
dynamical systems with many degrees of freedom, 
Josephson junction series arrays (JJSA) have been 


a subject of active research. After scaling the pa¬ 
rameters [Dominguez et al, 1991], the dynamical 
equations of an underdamped JJSA shunted by a 
resistive load, and subject to an rf-bias current 
I(t) = /(j c + 7 r fsin(ui r f£), [Hadley & Beasley, 1987; 
Hadley et al, 1988; Tsang et al., 1991; Tsang & 
Schwartz, 1992] are 

4>i + gf>i + sin <f>i+i L — idc + hi sin(fi r fr), 

( \ G V' i -1 AT W 

1 L = ctv(t) = — 2^ 99j, 1 = 1,..., N , 

V 3 =1 

where <f>i is the superconducting phase difference 
across the junction i. N is the total number of 
Josephson junctions or system size. Here, we use 
reduced units, with currents normalized by the crit¬ 
ical current, i = Ifl c \ time normalized by the 
plasma frequency u p t = r, with u> p = (2eI c /hC) 1 ^ 2 
and C the capacitance of the junctions; and volt¬ 
ages by rl c , with r the shunt resistance of the junc¬ 
tions. Il is the current flowing through the resis¬ 
tive load; g = ( h/2eCr 2 I c ) 1 / 2 = (3 C 1 ^ 2 , with j3 c the 
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McCumber parameter [McCumber, 1968; Stewart, 
1968]; v = Vtotal /N is the total voltage across the 
array per junction; a = rN/R, with R the resis¬ 
tance of the shunting load, represents the strength 
of the global coupling in the array; and the normal¬ 
ized rf frequency is fi r f = uj t f/u) p . Equation (1) ex¬ 
hibits rich spatiotemporal behavior, including phase 
locking, bifurcations, chaos, solitonic excitation, 
and pattern formation, breaking the law of large 
numbers and novel pseudo-Shapiro steps emerge in 
turbulence [Benz et al., 1990; Dominguez et al., 
1991; Dominguez & Cerdeira, 1995; Eikmans & 
van Himbergen, 1991; Free et al., 1990; Hebboul &; 
Garland, 1991; Kvale &; Hebboul, 1991; Lee et al., 
1992], However, to the best of our knowledge, the 
mechanism of the transitions among these dynam¬ 
ical phases, specially the transition from coherence 
to turbulence, has never been discussed. In this pa¬ 
per we study the interesting spatiotemporal inter- 
mittency, clustering bifurcation and clustering in¬ 
duced spatiotemporal chaos in the system (1). 

For a single cell (i.e. N — 1), the dynamical 
equation reduces to 

(j> + gcj) + sin 4> = z dc + z r fsin(f2 rf r), (2) 

with g = (1 + a)g. It is well known that Eq. (2) 
can exhibit chaotic behavior in the underdamped 
regime, i.e. g < 1 and fi r f < 1 [Ben-Jacob et al., 
1982; Bhagavatula et al., 1992; Huberman et al., 
1980; Iansiti et al., 1984; Jensen et al., 1984; Kautz 
& Monaco, 1985; Octavio & Raedi Nasser, 1984]. 
In Figs. 1(a) and 1(b) we show the bifurcation dia¬ 
grams, for g = 0.2, D r f = 0.8, as a function of i r f and 
id c respectively. In Fig. 1(a) with hie = 0.03, the fol¬ 
lowing points are to be remarked: First, the motion 
of the system (2) is period-1, then as i T f increases 
to a critical value 0.662, in the system takes place 
a period-doubling bifurcation to period-2. Second, 
as i r f continuously increases, the system undergoes 
a series of period-doubling bifurcation leading to 
a small scale region of chaos. At i T f « 0.832, 
this chaotic attractor suddenly expands, and is re¬ 
placed by a large scale chaotic motion. After the 
expanding transition the system acquires a rotat¬ 
ing motion, and the time-averaged voltage becomes 
nonzero. The bifurcation diagram as a function of 
Zdo with i v f = 0.61, is shown in Fig. 1(b). The bi¬ 
furcation behavior is essentially different from that 
of Fig. 1(a). As increases, the period-1 orbit 
first loses its stability, then a new period-2 solution 
arises via period-doubling bifurcation. The most 


(a) 




Fig. 1. Plots of <f> at t = nT(T = 2ir/Q, r f), with n being 
large enough to exclude the transient process. 

interesting and surprising point is that this period¬ 
doubling solution meets with an unstable period-2 
orbit (the dashed lines), and they suddenly disap¬ 
pear via inverse tangent (saddle-node) bifurcation 
as id c reaches a critical value ~ 0.035076. Be¬ 
yond this threshold, the behavior of the system is 
rotating and the motion is chaotic in a large scale 
region, and has the characteristic of type-I intermit- 
tency [Pomeau &; Manneville, 1980]. In Fig. 1(b), 
it is clear that another period-2 orbit appears via 
tangent bifurcation for id c near zero. Increasing 
idc, this period-2 solution first bifurcates into a 
small region of chaos through a series of continu¬ 
ous period-doubling bifurcations, then this chaotic 
motion coincides with the unstable period-2 orbit, 
and suddenly disappears due to a boundary crisis 
[Grebogi et al., 1983]. The two attractors form an 
interesting hysteresis phenomenon. In the following 
we investigate the complicated spatiotemporal dy¬ 
namics in JJSA and how it originates from that of 
a single Josephson junction. 

An important concept in a model for globally 
coupled systems is “clustering”. This means that 
even when the interaction between all elements is 
identical, the dynamics can break into different clus¬ 
ters, each of which consists of fully synchronized 








Rf-Driven Josephson Junction Series Arrays 1715 


elements. After the system falls in an attractor, we 
say that the elements i and j belong to the same 
cluster if </>, = <pj for all time. Therefore, the behav¬ 
ior of the whole system can be characterized by the 
number of clusters n c \, and the number of elements 
of each cluster (Mi, M 2 ,..., M„ cl ) [Dominguez k 
Cerdeira, 1995; Kaneko, 1989]. 

The simplest attractor of the system (1) is the 
spatially homogeneous configuration, so called co¬ 
herent state, i.e. 4>i (t) = 0(r), n c \ = 1, Mi = N. 
Linearizing Eq. (1) around the 4>{t) state 

N 

.. , (J _ , 

S<f>i + g5<fii + cos — g5<f>j = 0, 

3 =1 



i = 1,..., N. 

Introducing the following coordinates defined 
N 

3=1 

Y k = 5(f>k ~ 5<f>k+i, k = 1,..., N- 1 



irf 


After simple algebra, the critical stability bound¬ 
aries of this coherent state are determined by the 
following set of linearizing equations: 

X + gX + cos cf>X = 0, 

.. . (5) 

Y k + gYk + cos <f)Y k = 0, k = 1,..., N - 1. 


Fig. 2. (a) Bifurcation diagram for the homogeneous or co¬ 

herent state in a versus i t f plane with N = 128. In the shaded 
region the coherent state is unstable due to the clustering bi¬ 
furcation. (b) Bifurcation sequences plotted versus i r f for 
g = 0.2, fi rf = 0.8, idc = 0.03, <7 = 0.1 and N = 128. (c) The 
number of cluster, n c \ versus i T j for the state of Fig. 2(b). 


The first equation in (5) is nothing but the equa¬ 
tion obtained from linearizing the single cell case 
[Eq. (2)]. The second one characterizes the evolu¬ 
tion of the difference of two cell perturbations. The 
interesting point here is that the second one has 
the same structure as the other except that it has g 
instead of the renormalized g. Since the difference 
between the two is proportional to o, the system 
recovers the single cell scenario when a = 0. The 
critical boundaries of the coherent state in the a 
versus i r f parameter plane are shown in Fig. 2(a) 
with g = 0.2, i& c = 0.03 and N = 128. In the 
white region, the coherent state (motion in time 
may be regular and irregular) is locally stable, while 
in the shaded region, the coherent state loses its 
stability, and bifurcates to a multicluster state. As 
the coupling strength a decreases to zero, the in¬ 
stability regions collapse to the discrete bifurcation 
points for a single cell (a = 0). After the coherent 
state loses the stability, lots of multiclusters are cre¬ 
ated in the JJSA. A class of interesting states are 
multiclusters with a uniform distribution of junc¬ 


tions per cluster (i.e. M\ = • • • = M nd , with M i: 
being the number of elements in the ith cluster), 
and each cluster may have the same motion except 
for uniformly distributed phase shifts. We focus on 
this kind of states, a period-m state with k clus¬ 
ters will be called TmCk state, and N = k x n, 
n = 1, 2, 3,.... It often happens that m = k, then 
the dynamics of the TkCk state is reduced to 

k 

Ht) + g<j>(r) + sin </>(r) + | ^ 

— Idc “f ^rf sin(fi r fT) . (6) 

To investigate the clustering bifurcations in 
JJSA with nonzero coupling, we show the asymp¬ 
totic state of the system (1) in Fig. 2(b) as a func¬ 
tion of i T f with a = 0.1, N = 128 and the other 
parameters equal to those of Fig. 1(a). However, 
the bifurcation diagram is essentially different from 
that of Fig. 1(a). The T1C1 (coherent) state first 


g<i> 
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undergoes a cluster-period-doubling bifurcation at 
i r f s=s 0.624 to create a stable T2C2 state. By 
increasing i r f, the state undergoes further cluster- 
period-doubling bifurcations leading to spatiotem- 
poral chaos. Figure 2(b) is interesting due to the 
following novel features. First, we find a cluster¬ 
doubling induced period-doubling. The bifurcation 
point value is below the period-doubling condition 
for a single Josephson junction. Global coupling 
leads to cluster doubling at this parameter, which 
induces period doubling in time. Second, we find 
a cluster-doubling sequence 1-2-4 (and the induced 
period-doubling sequence). We expect that this 
clustering doubling cascade will proceed to a very 
large number of clusters and long periods. In our 
case this cascade is interrupted at k = 4 by a Hopf 
bifurcation, i.e. the modulus of another couple of 
complex eigenvalues is greater than one. Neverthe¬ 
less, the tendency of cluster doubling bifurcations 
leading to spatiotemporal chaos can still be seen 
in Fig. 2(c), where we plot number of clusters ver¬ 
sus i r f for the state described in Fig. 2(b). There¬ 
fore, we conclude that spatiotemporal chaos is made 
possible by clusterization, and call it “clustering 

(a) 

3 

2 

1 

"O- o 
-1 
-2 
-3 


0.0207 0.021 0.0213 



0.02 0.024 0.028 0.032 


idc 

Fig. 3. (a) Bifurcation sequences plotted versus -< c ic for g = 

0.2, firf = 0.8, j rf = 0.61, (7 = 0.6 and N = 128. (b) The 
critical boundaries among the T1C1 state, the T2C2 state 
and the turbulent phase in idc~cr plane with the parameters 
of (a). 


induced spatiotemporal chaos”. Moreover, these 
cluster-doubling sequences grow from the period¬ 
doubling sequences of the single cell due to the 
nonzero global coupling. As a decreases to zero, 
the clustering-doubling sequences is identified as the 
period-doubling sequence of the single cell. If the 
period-doubling sequence of the single cell is bro¬ 
ken off, then the character of the clustering bifur¬ 
cation in JJSA also changes suddenly. This can be 
clearly seen in Fig. 3(a) which shows the asymp¬ 
totic state of the system (1) along the axis, with 
a = 0.6 and the other parameters are the same as 
those of Fig. 1(b). The T1C1 state first undergoes a 
cluster-period-doubling bifurcation at ~ 0.02087 
to create a stable T2C2 state. However, since the 
period-doubling period-2 solution in the single cell 
[see Fig. 1(b)] is destroyed by the inverse saddle- 
node bifurcation by increasing id c , the T2C2 state 
in the JJSA is suddenly destroyed by the spatiotem¬ 
poral intermittency transition near tdc = 0.02143. 
In Figs. 2(b) and 3(a), first we run Eqs. (1) to get 
the coherent (period-1) state from random initial 
conditions, then we compute Eqs. (1) by gradually 
increasing the parameter value (ij c or i r f) and by 



-1.3 


-i-5! i uu in u uiuumuuu! uuu iu 

- 1.6 - 

0 20 40 60 80 100 



i 


Fig. 4. Snapshot of the asymptotic solution of the sys¬ 
tem (1) after the transient process for g = 0.2, ?' r f = 0.61, 
firf = 0.8, a = 0.6, and N = 100. (a) idc = 0.0206, (b) and 
(c) are two successive snapshots for i& c = 0.0210. 
















using the final state for the previous parameter 
value as the initial state for the new parameter 
value, in this way we can surely get clusters with a 
uniform distribution of cells for all cluster-doubling 
cascades. Figure 3(b) shows the phase diagram 
among T1C1 state, T2C2 state and the turbulent 
phase in the a versus i& c plane. The two criti¬ 
cal transition curves in Fig. 3(b) are obtained by 
the numerical simulation of the system (1). It is 
clear that the regime of the T2C2 state is very nar¬ 
row. Figures 4 shows the snapshots of <f> for the 
T1C1 state and T2C2 state after a long transient 
process. The features of coherence and two-cluster 
are clearly observed in Figs. 4(a), and 4(b)-4(c), 
respectively. The most interesting phenomenon is 
that the system suddenly evolves to a very compli¬ 
cated rotating motion as i& c is increased beyond a 
critical value (idc ~ 0.02143 for a = 0.6), i.e. af¬ 
ter the T2C2 state loses its stability. The system 
falls in a large n c \ ~ N clusters motion with all 
Mj small. Figure 5 shows the space-time evolu¬ 
tion after a very long and complicated transient 
process for = 0.0215 and a = 0.6. The tur¬ 
bulent character of the motion is very clear. The 
evolution of (f>i (the first junction) is displayed in 
Figs. 6(a) at the same parameters values as those of 
Fig. 5. The motion displays periodic behavior (2P) 
for a long time, it is suddenly interrupted by large 



Fig. 5. The time-space evolution of the system (1) for 
idc = 0.0215 and the other parameters the same as those of 
Fig. 4. The plots are at t = nT after some transient process, 
where T is the same as in Fig. 1. The features of turbulence 
are clearly observed. 
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Fig. 6. The evolution of <f>\ and <pi — <j >2 with the same pa¬ 
rameters as those of Fig. 5. The features of spatiotemporal 
intermittency are clear. 

bursts and quickly resumes the periodic fashion. 
The similar features of the difference (pi—4>2 are also 
displayed in Fig. 6(b). As i,\ c is far from the criti¬ 
cal value, more and more random large bursts take 
place more frequently. Although this behavior is 
similar to the characteristic of well-known intermit¬ 
tency, which were investigated in low-dimensional 
systems [Pomeau & Manneville, 1980], it is an es¬ 
sential type of spatiotemporal intermittency, which 
has not been found before in the rf-driven JJSA 
or other high dimensional globally chaotic systems. 
The above features do not depend on the specific 
cell and the number of cells. The spatial variable 
(N > 2), the dynamics of a single cell and the global 
coupling are of crucial importance for this interest¬ 
ing phenomenon. 

In conclusion we analyzed the complex spa¬ 
tiotemporal dynamics of the rf-driven JJSA. Clus¬ 
tering bifurcation, clustering induced spatiotem¬ 
poral chaos and spatiotemporal intermittency are 
shown to appear in these systems. The spa¬ 
tial variable, the dynamics of a single cell and 
the global coupling are of crucial importance for 
the existence of these interesting spatiotemporal 
phenomena. 
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In this paper, we are devoted to the problem of escaping from a potential well which is present 
in a great number of physical situations. We use the Helmholtz oscillator as a model for those 
situations and consider the behavior of the oscillator under an additional driven perturbation. 
The Melnikov analysis reveals it as an adequate method. Some comparisons are made with the 
perturbations of the oscillator on the linear and quadratic terms. 


1. Introduction 

The problem of escaping from a potential well is 
present in a great number of different physical sit¬ 
uations: The orbits described by a photon near a 
Schwarzschild black hole [Misner et al., 1973], the 
capsizing of a boat under trains of regular waves 
[Thompson, 1989] and others (see [Chacon et al., 
1996]). 

In this paper we continue developing the idea 
first introduced in [Chacon et al., 1996] where the 
Helmholtz oscillator was considered as a model to 
treat most of the above problems. 

In all those situations a common characteristic 
can be observed. Before the escape from a potential 
well, some chaotic transients of an unpredictable 
length can be produced and they can be observed 
in the orbits which start in the chaotic regions of the 
phase space, one of those being the region closed to 
separatrices [Chacon et al., 1995]. 


In [Chacon et al., 1996] we applied to the 
Helmholtz oscillator the technique of weak paramet¬ 
ric modulations to the quadratic and linear terms 
in the equation of the oscillator. The technique 
attains good results in the task of the inhibition 
of the chaos which was present in the unperturbed 
oscillator. 

The Helmholtz oscillator was chosen since it is 
a good example and a representative model of the 
behavior of the situations mentioned above. 

In its integrable form its equation can be given 
by: 

x"{t) — x(t ) + (3x(t ) 2 = 0 

where x denotes the displacement and (3 > 0. 

When we introduce a weak perturbation in the 
integrable expression, the equation takes the form 

x"(t) — x(t) + (3x(t) 2 = — 5x'(t) + 7 sin(u;£) 


‘This work has been supported by Spanish D.G.I.C.Y.T. PB95-1004 grant and also by COM-20/96 MAT from Direction 
General de Universidades, Comunidad Autonoma de Murcia. 
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with 5 > 0 and 7 <C 1 . c 0 , 5 and 7 denote the nor¬ 
malized frequency, damping coefficient and driving 
term amplitude. 

In this paper we are devoted to study the be¬ 
havior of the oscillator under an additional driven 
perturbation given by: 

x"(t ) — x(t) + (3x(t) 2 = —8x'(t) + 7 sin(cjf) 

+ cry sin(fl£ + tp) ( 1 ) 

where 0 < a <C 1 and with O, and p we denote the 
normalized frequency and the initial phase of the 
perturbation. 

We will compare the results with those obtained 
in [Chacon et al ., 1996] where perturbations of the 
quadratic and linear terms were considered. As a 
conclusion we will have that in some cases it is more 
advantageous to use the additional driven term (low 
frequencies) while in others the use of perturbations 
in the quadratic and linear terms (high frequencies) 
is better. 


tables [Gradshteyn et al., 1994] we obtain: 

M(to) = — C — A cos(u>to) — B cos(fl£o + <p) 
where 

66 67T7 o . , . 

C = 5^2 A = -/T uj cosh(mj) 

B = ft 2 cosh( 7 ro>) 

Note that A, B and C are positive numbers in all 
the interval of the parameters. 

It is well known that the Melnikov’s function 
M(to) gives a measure of the distance between 
the perturbed stable and unstable manifolds in the 
Poincare section at a particular time to■ When 
M(t 0 ) has a simple zero, then we have a homo¬ 
clinic bifurcation and therefore we have the possi¬ 
bility of the appearance of chaos. As a consequence 
M(t 0 ) = 0 is a necessary condition for such an ap¬ 
pearance. But if M(t 0 ) has never a zero then this is 
a sufficient condition for the opposite effect of the 
inhibition of even the transient chaos. 


2. Melnikov Analysis 

When detecting chaos in a dynamical system, 
one analytic procedure of great interest is the 
Melnikov’s method [Wiggins, 1990], It gives a 
useful criterion for detecting the presence of ho¬ 
moclinic or heteroclinic orbits. Nevertheless this 
procedure is limited in two senses: It is an ap¬ 
proximative method of first order and is only valid 
for orbits which start in points sufficiently close to 
separatrices. 

In the situation under consideration it gives 
good results. If we consider the phase space cor¬ 
responding to Eq. (1), the separatrix of the system 
without perturbation has the following equations: 

xo(t) = 3/2/3 sec h 2 (t/ 2) 

x' 0 {t) — —3/2/3 sec h 2 (t/2)tgh(t/2) 


The Melnikov’s function associated to such an 
orbit is: 

/ OO 

Xq ( t)dt 

-OO 


+ 7 


/ OO 

x' 0 (t) sin[w(f + to)]dt 

-OO 


3. Inhibition of Chaos 


When a = 0 then M(fo) = — C — A cos(uAo) and we 
have a situation of chaotic escape if A — C = <3 > 0 
since in this case the Melnikov’s function has simple 
zeros and it is equivalent to the previous condition 
d> 0 . 

If a 7 ^ 0 then B > d, that is A — B — C < 0. 
This last relationship is a necessary condition for 
M(to) to have the same sign, in this case M(to ) < 0 
for every to and can be written as 


a > 


,2 r 


sinh(7rfi) 


C\ uj 

1 - 7 it where R= 777 , . . , . , 

AJ II 2 Lsmh(7Rj)J 


For general values of and 0 < < 2 tt 

the above condition is not sufficient to assure that 
M(to ) < 0 for every to- Now we can state the 
main result on inhibition of chaos for the Helmholtz 
oscillator. 


Theorem 1. Let LI = puo with p a positive in¬ 
teger and such that p = —^ ~n+i^ satisfied for 

some positive integers m and n. Then M(fo) always 
has the same sign (M(fo) < 0) for every to if and 
only if 


+ cry 



x'q (t) sin[fi(f + to) + <p\dt 


After some computations and using the integral 


where 


C*inin<C¥ Ct max 


C*min — 



CkmRx 


R 
P 2 
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One observation on the theorem is that the 
length of the interval [a m j n , a m ax] could be small, 
in fact this is the case if we observe the val¬ 
ues of the parameters for which we have chaotic 
escape. 

To prove the theorem we need some technical 
lemmas. 

Lemma 2 . Let pu = qCl for some positive in¬ 
tegers p and q. Then there exists satisfying 
— cos(cufo) = cos(flfo + <p) = 1 if and only if 
q = f or some integers m and n. 

Proof. If — cos utQ = 1 then there is an integer n 
such that 

ct>fo = ( 2 n + 1 ) 7 r ( 2 ) 


A sin 


u>to + 



B sin(fito + o') 


= A — B cos 



> A-B 


since Cl and cu are incommensurable and 
cos(2tt£) < 1. 

In an equivalent way we will take a = | — 
37 t^ = § (1 — 17 ) when d < Cl. Let a 7 ^ | — 37 r^ = 
|(1 — ^). Then taking ?o = 77 we obtain 

A sin futo + — B sin(fI?o + o) 


= A — B sin 




> A-B 


since in this case a 7 ^ 5 — 37 t —. ■ 

1 Z U) 


On the other hand if cos(f2fo + p) = 1 then there 
is m 6 Z such that CUq +.tp = 2mrc, that is, 
f Uq = 2 m 7 r — ip. Using (2) we obtain 

Cl _ 2rmr <p p 

u (2n + l) 7 r ( 2 n+l) 7 T q 

Therefore 

p 2mir — p 2m — % 
q (2 n + l)ir 2n + 1 

Lemma 3. Let Cl and u be incommensurable, that 
is, 77 is an irrational number. Then there exists tn 
such that 

—B cos(Clio + <p) — A cos(u>to ) > A — B 

Proof. Using the change of variable <p = a + ^ we 
have 

—B cos(flto + <p) — A cos (cut 0 ) 

= A sin (uto + ^p\ - B sin(J2fo + cr) 


Lemma 4. Let g{r\ p, q) = where r € 

R and p and q are positive integers. Then g is fi¬ 
nite if and only if q = 1. Moreover, we have that 
0 < g(r\ p, 1) < p 2 , where r 6 (— 00 , 00 ). 

Proof. First we will compute lim T _ > 2 7 r/ <?( r ; P> Q) 
where l G Z. A necessary condition to get g 
bounded is that its numerator has zeros in the same 
points that its denominator. 

Consider first the case 5 / 1. The zeros of 
the numerator and denominator of g are r = ~2ln 
t = 27rs where I £ Z and s£Z. It is evident that 
with arbitrary positive integers p, q(q 7 ^ 1 ) it is not 
always possible to find an integer l such that s = 1^ 
where s is an arbitrary integer. 

The above limit attains the value p 2 and since 
2 

, using finite induction over 
p we immediately obtain 0 < g{r\ p, 1) < p 2 . 

Proof of Theorem 1. We ask under which condi¬ 
tions the condition a > (l\— ^)R is also a suffi¬ 
cient condition to obtain M(fo) for every to- It is 
immediate to observe that a sufficient condition is 
given by 

A — B > —B cos(ff£o + <p) — A cos(cufo) (3) 


g{r\ p, 1 ) = 


sin(pr/2) 

sin(r/2) 


We are looking for a to for which 

A sin ^u7?o + - 77 ^ — B sin(fI?o + cr) > A — B 

Let a = | — 37 Tj = |(1 — if u > Cl. Then 
taking to — 7 we obtain 


Now we are looking for values of cu, U and p in 
order that Eq. (3) be satisfied for every to- Lemma 3 
says that if Cl and cu are incommensurable, this con¬ 
dition is not satisfied. Therefore a resonance condi¬ 
tion pw = qCl is necessary with p, q positive integers. 
In such case Lemma 2 supplies a sufficient condition 
for Eq. (3) to be held in a infinite number of values 
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0.00 0.40 0.80 1.20 


00 

Fig. 1. Function Aq(u) = a m ax — a m in versus u for <5 = 
const., (3 — const., 7 = const, and 11 == pio: (a) p = 1; 
(b) p = 2; (c) p = 3. 

of to. In such a case (3) is written as 

A 1 — cos(flfo + <p) 

B ~ 1 + cos(wfo) 

1 — cos ^-[uto — (2n + 1 )tt] + 2 m7r^ 

1 + cos(<ufo) 

Then 


Act achieves a minimum value as a function of the 
frequency. 

4. Comparisons with the Quadratic 
and Linear Parametric 
Perturbations Cases 

In [Chacon et al., 1996] the Helmholtz oscillator is 
perturbed by weak parametric modulations of the 
quadratic (5a) and linear (5b) terms in the following 
way: 

x"(t) = x(t) — P[1 + r} sin(ffif + pi)\x(t) 2 

— 5x'(t ) + 7 sin(u;f) (5a) 

x"(t) = x(f)[l + rf sin (fl't + ip')\ 

— f3x(t) 2 — Sx'(t) + 7 sin(wt) (5b) 

obtaining similar results as above. Nevertheless it 
is of great interest to compare those results from 
the point of view of inhibition of chaos with the sit¬ 
uation considered in this paper. The results of this 
comparisons are shown in Figs. 2 and 3. In both 
pictures we represent the quotient ^ Aa( ^) versus 
uj for several values of p. In Fig. 2 we have 
A 7](lo) _ IO 7 (3 

Aa(w) (p 2 u 2 + 1 )(p 2 uj 2 + 4) 
and in Fig. 3 

Ap'(uj) _ 2q f3 

Aa(u) p 2 u 2 + 1 


1 P T 

1 — cos - 


B 


> 


1 — cos r 


But since A 


R 

— > 
a 


1 — cos 


w 


pr 

Q 


cos r 


with t - coto — (2 n + l) 7 r 

( 4 ) 

= § we obtain 

L sinh(7ro;;J a 


with t = uto — (2 n + 1)tt 


Finally if q = 1, Lemma 4 gives a condition for 
(4) to be satisfied for every r (that is for every to)- 

' ]• 


This condition is a < where R — ^ 


sinh(7rn) ] 
sinh(7ro;) 


Remark 1. Suppose we have a set of parameters 
that satisfy the hypotheses of Theorem 1. As p in¬ 
creases, the interval of escape (a m i n , a m ax) shrinks 
rapidly. Figure 1 shows a plot of the width of the 
inhibition interval as a function of uj. For each p. 



Fig. 2. The function 77 versus w f° r several values of 

p: (a) p = 1; (b) p = 2; (c) p = 3. 







(1 /y P)(Ar( , (co)/Aa(co)) 



Fig. 3. The function Aa(<^) versus u for several values of 
p: (a) p = 1; (b) p = 2; (c) p = 3. 

From these two pictures we can see that for a 
fixed value of p, an additional damped perturbation 
of the oscillator is more efficient (in the frequency 
domain) than perturbations of the quadratic and 
linear terms for low frequencies u for the inhibi¬ 
tion of chaos. The opposite effect is obtained for 
large frequencies u. It can be observed also that for 
oj = p we obtain Atj = Aa and At/' = Aa. 

5. Conclusions 

The application of an additional driven pertur¬ 
bation to the Helmholtz oscillator is an efficient 
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method to suppress the chaotic escape from a po¬ 
tential well. We have estimated the range of pa¬ 
rameters where the inhibition of chaos is possible 
through a systematic application of the Melnikov’s 
method. 

When we compare the results with those ob¬ 
tained by parametric perturbations of the quadratic 
and linear terms we obtain that for a fixed value of 
the resonance parameter p, at low frequencies uj the 
effectiveness of the suppression of chaos by an addi¬ 
tional driven term is greater than in the other two 
methods. At large values of the frequency we obtain 
the opposite situation. 
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The mathematical modeling of biological systems has proven to be a valuable tool by allowing 
experiments which would otherwise be unfeasible in a real situation. In this work we pro¬ 
pose a system of nonlinear differential equations describing the macroscopic behavior of the 
cardiac conduction system. The model describes the interaction between the SinoAtrial and 
AtrioVentricular node. Its very simple structure consists of two nonlinear oscillators resistively 
coupled. 

The numerical analysis detects different kinds of bifurcations whose pathophysiological 
meanings are discussed. Moreover, the model is able to classify different pathologies, such 
as several classes of arrhythmic events, as well as to suggest hypothesis on the mechanisms 
that induce them. These results also show that the mechanisms generating the heartbeat 
obey complex laws. The model provides a quite complete description of different pathological 
phenomena and its simplicity can be exploited for further studies on the control of cardiac 
dynamics. 


1. The model 


We propose to describe the interaction between the SinoAtrial node and the AtrioVentricular node by 
modeling them as two-coupled nonlinear oscillators. The model equations can be associated to an equivalent 
electronic circuit, depicted in Fig. 1. 

The system of differential equations describing the model is then: 


x\ = 77 - x 2 
w 

x 2 = —— [x\ + g(x 2 ) + R(x 2 + x 4 )] + A cos(27r ft) 
L 1 

1 

X3 X4 

X4 = - 7 - [xz + /(x 4 ) + R(x 2 + z 4 )] 

l 2 


( 1 ) 
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Fig. 1. Equivalent electrical circuit describing the model 
constituted by two oscillators: SA and AV. The parame¬ 
ters have the following values: R = 0.9 ft, C\ = 0, 167 F, 
Li = 0, 033 H, C 2 = 0, 45 F, L 2 = 0, 018 H. 


where 


f{x 4 ) = -X 4 + ^£4 


with 


h(x 2 ) = < 


1 3 

-£ 2 + -I 2 

+ h(x 2 ) 

2 1 
~ x 2~4 

Nl < \ 


1 

-X 2 



1 

X 2 

X2 <2 


( 2 ) 

( 3 ) 


( 4 ) 


We assume that x 2 and £ 4 , describe, respec¬ 
tively, the action potentials of the SA and AV nodes 
[di Bernardo et al ., 1998]. 

In the past ten years, even more if we consider 
the pioneering work of van der Pol [van der Pol 
& van der Mark, 1928], many attempts have been 
made to describe the macroscopic heart activity by 
means of mathematical modeling. All of these mod¬ 
els can be classified in two categories: ( 1 ) discrete 
time models based on circle maps [Guevara & Glass, 
1982; Bub & Glass, 1994; Honerkamp, 1983], and 

( 2 ) continuous time models based on limit cycle os¬ 
cillators [van der Pol & van der Mark, 1928; West 
et al., 1985]. The time continuous model presented 
in this paper differs somewhat from others of its 
class in four main features: 


(1) Each oscillator is described by a two- 
dimensional system, so that the system of or¬ 
dinary differential equations [Eqs. (1)] describ¬ 
ing the complete model is 4(2 + 2)-dimensional, 
thus enabling complex behaviors such as ape¬ 
riodic and chaotic solutions. In the model of 
West [West et al., 1985], on the contrary, the 
oscillators are monodimensional, consequently, 
the complete model equations are 2(1 + 1 )- 
dimensional. 

(2) Each oscillator is able to oscillate on its own, 
that is, without an external forcing, thus re¬ 
specting the autoexcitatory nature of both SA 
and AV node cells (whereas, in the West model, 
a voltage generator is needed). 

(3) The coupling between the two oscillators is bidi¬ 
rectional, therefore different from the van der 
Pol model [van der Pol & van der Mark, 1928], 
where it is monodirectional. 

(4) Unlike other models of the same kind, the two 
oscillators are not identical as they should take 
into account the different physiological behavior 
of SA and AV nodes. 

In short, our model put together the features of 
the van der Pol model with the West model one in 
a compact and simple set of equations. 


1.1. The choice of parameters values 

Parameters values were fixed respecting the follow¬ 
ing physiological constraints: 

(1) The shape of the model outputs x 2 (SA node 
Action Potential) and £4 (AV node action po¬ 
tential). They depend on the ratios £\ = C\/L\ 
and e 2 = C 2 /L 2 . Figure 2 shows an example of 
this behavior. A change in the value of C\ corre¬ 
sponds to different slopes of the rising wavefront 
of £ 2 . This is similar to the behavior experimen¬ 
tally observed in SA node cells as their depo¬ 
larization frequency varies [DiFrancesco, 1995]. 
This change of the slope is the mechanism used 
by the Central Nervous System to control the 
heart rate and is achieved by changes in ionic 
currents and in cellular membranes permeabil¬ 
ity [DiFrancesco, 1995]. 

(2) We chose the physiological shooting value of 
the SA node. Indeed, the AV node is also a 
pacemaker site: in pathologies that prevent the 
SA node from depolarizing (e.g. Sinus Arrest 
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Fig. 2. Waveforms of X 2 for two different values of C'i. A 
larger value of C\ corresponds to a smaller intrinsic frequency 
of the SA node oscillator. 




Fig. 3. Length of the simulated P-R intervals (top diagram) 
and the period of the first oscillator, that represent the SA 
node (bottom diagram), as a function of the coupling resis¬ 
tance. The parameters are chosen as in Eq. 2. 


event), the AV node becomes the dominant 
pacemaker of the heart. However its shoot¬ 
ing frequency is lower than the SA node one 
so that the heart rate also slows down. This 
phenomenon is known as junction rhythm. It 
has been found in experiments on canine hearts 
that the ratio between the physiological SA and 
AV nodes shooting frequencies can be consid¬ 
ered equal to Tsa/Tav = 2/3 [West et al., 
1985], 

(3) In the normal ECG signal, the time delay be¬ 
tween the P wave (onset of atrial contraction) 
and the R wave (onset of ventricular contrac¬ 
tion) is in the range [0.12 s, 0.25 s]. In our 
model the P-R interval is the delay between 
the maximum value of x 2 (SA Action Poten¬ 
tial) and the maximum value of x 4 (AV Ac¬ 
tion Potential). This time interval depends 
mainly on the coupling resistance R. Figure 3 
shows both the obtained P-R and Tsa val¬ 
ues as R varies. We observe an acceptable 
range of P-R intervals for R > 0.8 ft. We set 
R = 0.9 f2 so that the corresponding value of 
Tsa is physiological and approximately equal 
to 80 b/m. 

The parameters values satisfying the above 
mentioned properties are: 

Ci = 0.167 F Li = 0.033 H C 2 = 0.45 F 

(5) 

L 2 = 0.018 H R = 0.9 n 


C'i, Li 

drive the frequency of the SA node 

C 2 , l 2 

drive the frequency of the AV node 

R 

drives the coupling between the 


two oscillators, as well as, the 


P-R time delay 


The absolute refractory period of the first oscil¬ 
lator (representing the SA node) with the parame¬ 
ters as in Eq. (5) is approximately 0.25 s, while for 
the second oscillator (representing the AV node) is 
approximately 0.15 s. This refractory periods where 
calculated by forcing the two oscillators, in uncou¬ 
pled condition, with a square impulse with an am¬ 
plitude of 5A and a duration of 10 ms. 

Technical details together with the extensive 
model description can be found in [di Bernardo 
et al., 1998]. 

2. Equilibrium Point Bifurcations 

With the parameters as in Eq. (5) the system 
has got only one unstable equilibrium point at 
(aq, x 2 , X 3 , £ 4 ) = (1/4, 0, 0, 0), where we set x\ — 
Vi, x 2 =Ii,x 3 = V 2 , x 4 = I 2 . 

We want to study the behavior of our model 
when the coupling resistance R and the capacitance 
Ci vary. 

The R parameter models the coupling 
“strength” between SA and the AV nodes. When 
R increases the coupling between the two oscilla¬ 
tors increases, because a smaller current will flow 
through the resistance. Changing the value of 
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Fig. 4. Plot of the real parts of the Jacobian matrix eigenvalues, evaluated at the equilibrium point, when C\ varies for 
different values of the coupling resistance R. When [0 < R < 0, 5] all four eigenvalues have positive real parts (A). When 
[0, 5 < R < 1] the real parts of two eigenvalues change their sign, and for small values of C\, they are all 4 greater than 0 (B). 
When [/? = 1] no eigenvalues change sign (C). For small values of C\ all four eigenvalues’ real parts become equal to zero. 
When [it > 1], the real parts of two eigenvalues change their sign, and for small values of C\ become all 4 less than 0 (D). 


C\ means to increase or decrease the depolar¬ 
ization frequency of the SinoAtrial node. The 
SA frequency increases when C\ decreases and 
vice-versa. 

To study the possible bifurcations of the equi¬ 
librium point of our system as R and C\ vary, we 
calculate the eigenvalues of the Jacobian matrix of 
our system at the equilibrium point. Figure 4 sum¬ 
marizes the results of our analysis. When R < 1, 
the equilibrium point remains always unstable, thus 
no interesting behavior is found. However when 
R = 1, there exist a value of C\ for which all four 
eigenvalues become pure complex numbers. This 
suggests that a double Hopf bifurcation takes place. 
We will further show this bifurcation by using nu¬ 
merical simulations. When R > 1, two eigenvalues 


become pure complex numbers, for a specific value 
of Ci. Thus we can expect the occurrence of a Hopf 
bifurcation [Kuznetsov, 1995]. 

A numerical simulation confirms what we hy¬ 
pothesized above. Figure 5 shows the results for 
double Hopf bifurcation. By decreasing R from 1.1 it 
to 0.9 Q with Ci = 0.0027 F, the stable fixed points 
become unstable and a torus is generated. 

Figure 6, instead, shows the Hopf bifurcation. 
We notice that this is a subcritical Hopf bifurcation 
for which the unstable limit cycle folds back and 
becomes stable. 

This kind of scenario (Hopf bifurcation plus fold 
bifurcation) is typical for biological systems [Glass 
k Mackey, 1988; Seydel, 1994] and is known as hard 
loss of stability of the equilibrium point. 
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(a) 



(b) 

Fig. 5. (a) Projection of the four-dimensional trajectory falling in the equilibrium point on the (X 3 , * 2 , X 4 ) phase space. 

R = 1.1 0, Ci = 0.0027 F, and the other parameters as in Eq. (2). (b) Projection of the torus in the (X 3 , X 2 , * 4 ) phase space. 
R — 0.9 Q, Ci = 0.0027 F, and the other parameters as in Eq. (2). 



A physiologic interpretation of these bifurcations 
is possible: when the intrinsic frequency of the 
SA node increases (that is C\ decreases) beyond 
a critical value, the AV node may not be able to 
shoot at the high frequency dictated by the SA 
node anymore (this corresponds, on the bifurca¬ 
tion diagram in Fig. 6, to the bistability region). 
The greater the increase of the SA node frequency, 
the higher the likelihood that the heart would stop: 
the bistability region ends and the only remaining 
stable attractor is the critical point. 

3. Simulation of Cardiac Arrhythmias 

In this section we will study the genesis of 
some cardiac arrhythmias that occur when some 


Fig. 6 . Subcritical Hopf bifurcation. The bistability region 
is shown as C 1 varies. The minima and maxima of X 2 are 
shown for different initial conditions, after a time interval 
long enough for the transient to be considered over. The 
unstable limit cycle is not shown. R = 1.1 fi. 
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parameters of the system change. We will also force 
the system on the SA node oscillator using a sinu¬ 
soidal voltage generator. The external generator 
action should model the influence of ectopic pace¬ 
makers in the atria that disturb the activity of the 
dominant pace-maker of the heart (the SA node). 

3.1. Unforced system 

The analysis for the unforced system begins by 
varying the coupling resistance R while the other 
parameters are constant with values as in Eq. (5). 

We are interested in which kind of phase¬ 
locking between the two oscillators occurs. As a 
matter of fact, different phase-locking behaviors 
correspond to different kind of arrhythmias. In 
order to classify the solutions, a specific rotation 
number p is defined [di Bernardo et al., 1998]. In¬ 
tuitively p represents the average number of times 
the SA node depolarizes for a single depolarization 
of the AV node. For example p = 1.5 means that, 
on average, for every three depolarizations of the SA 
node, two depolarizations of the AV node occur. We 
calculated the rotation number p, for different val¬ 
ues of the coupling resistance R. For R > 0.11 tt the 
two oscillators are 1:1 phase-locked. Every time SA 
depolarizes, there is one AV node depolarization. 
As R decreases a series of subharmonic bifurcations 
appear. Some of these solutions are aperiodic, but 
all of these resemble just one type of arrhythmia 
known as 2° AV block of the Wenckebach type. 

A similar analysis is performed when R is fixed 
and C\ varies. This corresponds to varying the fre¬ 
quency of the first oscillator (SA node). It hap¬ 
pens that the rotation number assumes only integer 
values n, meaning that the phase-locking between 
the two oscillators is of the n: 1 kind. These so¬ 
lutions resemble arrhythmic episodes known as n:l 
AV blocks. Thus our model is able to simulate dif¬ 
ferent kinds of arrhythmia depending on which pa¬ 
rameter, R or C\ , is changed. However we are also 
interested in the study of solutions that could re¬ 
semble atrial and ventricular fibrillation. There¬ 
fore we searched for chaotic solutions in the R- 
Ci parameter space. We found a region of the 
parameter space (0.967 Vl < R < 1 and 
0.012 F < C\ < 0.015 F) in which a series of 
period doubling bifurcations, leading to a chaotic 
attractor, takes place. From the pathophysiological 


point of view, this result is important because it 
shows that for high depolarization frequency of 
the SA node [C\ = 0.012 F, period of SA oscil¬ 
lator ^=300 beats per minute), a chaotic behavior, 
such as fibrillation can occur. 

3.2. Forced system 

We set the parameters as in Eq. (5), but R = 
0, 11 This is done as more complex dynamics 
occur at the edge of the 1:1 phase-locking region, 
found for this value of R. 

In order to identify different kind of arrhyth¬ 
mias, we classify the solutions by calculating 
Poincare maps. The solution is sampled at inter¬ 
vals equal to the forcing function period. Using 
these maps it is possible to define a rotation num¬ 
ber pi which has the same meaning given in the 
previous paragraph. Moreover, it is possible to cal¬ 
culate, when the solution is periodic, the number of 
periods n and m of the first and second oscillator, 
contained in one period of the forcing function. 

The solutions obtained as the frequency and 
the amplitude of the forcing function vary, resem¬ 
ble both Wenckebach rhythms and AV blocks and 
also atrial bigeminy like and chaotic solutions. The 
rotation number of Wenckebach-like solutions as¬ 
sumes rational values; AV blocks like solutions, in¬ 
stead, have integer rotation numbers. The solu¬ 
tions that seem to be chaotic are generated at high 
frequency and low amplitude values of the forcing 
function. This result suggests that an ectopic pace¬ 
maker (modeled by the forcing function) can initiate 
a chaotic event in the heart muscle. 

Atrial bigeminy episodes were obtained for 
small values of the frequency of the forcing function, 
smaller than half of the free shooting frequency of 
the SA oscillator. 

By increasing the coupling “strength” (R from 
0.11 to 0.13 fl), the kind of arrhythmia simu¬ 
lated by the solution (atrial bigeminy) remains the 
same. Only the R-R interval series change. 1 For 
R — 0.13 ft the R-R series, shown in Fig. 7, seems 
to be modulated by an aperiodic function, whereas 
for R = 0.11 ft, R-R interval series remains periodic 
and alternates between two fixed values. 

This behavior is significant as it has been ex¬ 
perimentally found [Babloyantz & Destexhe, 1988] 
that the normal heart is not a perfect periodic 


’The time interval between two consecutive QRS complexes is known as R-R interval. We calculated them in our model as 
the time intervals between two consecutive maxima in x±. 









Fig. 7. Intervals between two successive X 4 maxima during 
atrial bigeminy event. Forcing function amplitude A = 1.1 V, 
frequency / = 0.4 Hz and coupling resistance R = 0.13 fl. 
The other parameters as in Eq. (5). 

oscillator. Its dynamics, indeed, are very rich and 
complex. In particular R-R interval series have 
been extensively investigated suggesting the pres¬ 
ence of nonlinear mechanisms of chaotic nature con¬ 
tributing to their generation and control [Guzzetti 
et al., 1996]. 

4. Conclusions 

The simple nonlinear model (someone can say raw!) 
proposed in this paper allows the characterization 
of several pathophysiological conditions that could 
occur or affect the real heart. 

It is able to reproduce many episodes like ar¬ 
rhythmias or other alterations in the normal heart 
rhythm. This is for example the case of cardiac 
arrest event. It can be explained through a transi¬ 
tion marked by a Hopf plus a Fold bifurcation when 
the frequency of the driving oscillator reaches values 
that are too high. 

Moreover in many regions of the parameters 
space, bifurcation diagram allows the interpretation 
of heart dynamics in terms of generating mecha¬ 
nisms. In this way, the analysis of the model by 
means of the powerful tool introduced by the bifur¬ 
cation theory applied to biological systems not only 
produces a new description of the heartbeat gener¬ 
ation but also helps in the identification of possible 
causes for these events. 

The results of the analysis could lead to ap¬ 
plications in heart rate prediction and control by 
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exploiting one of the features of the model such as 
the reduced number of parameters and their strong 
correspondence to the physiological cardiovascular 
system behavior. 
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Rings of chaotic oscillators coupled unidirectionally through driving are studied. While synchro¬ 
nization is observed for small sizes of the ring, beyond a certain critical size a desynchronizing 
transition occurs. In the two examples studied here the system exhibits a transition to periodic 
rotating waves for rings of Lorenz systems, while one finds a sort of chaotic rotating waves 
when Chua’s circuit is used. 


1. Introduction 

Synchronization phenomena are pervasive in na¬ 
ture [Winfree, 1980; Strogatz & Stewart, 1993] and, 
thus, many studies have been carried out, focusing 
particularly on limit-cycle oscillators. Less intuitive 
is probably the finding that chaotic systems may 
be made to get in synchrony [Fujisaka & Yamada, 
1983; Pecora & Carroll, 1990], as chaos has been 
described as a situation in which a system gets out 
of synchronization with itself [Tang et ai, 1982]. In 
the present work we shall use the synchronization 
method introduced in [Giiemez &; Matfas, 1995], 
that amounts to a generalization of the method in¬ 
troduced by Pecora and Carroll (PC). The idea is 
to avoid partitioning the response system in sub¬ 
systems, introducing, instead, the driving signal at 
a particular place of the response system, i.e. with¬ 
out reducing the size of the latter. This property 
is particularly useful in the case that we are inter¬ 
ested here: the design of arrays of coupled chaotic 
oscillators, as all the units in the array will be of the 
same type (will have the same dimension), without 


reducing the richness in possible dynamical behav¬ 
iors of the system. This method has been used 
before and synchronization waves in linear arrays 
of chaotic oscillators have been obtained [Sanchez 
et al., 1997]. 

In the present work, we shall consider rings of 
coupled chaotic oscillators. These geometries may 
be relevant in a biological context, like in morpho¬ 
genesis [Turing, 1952] or in the context of neural 
systems. Thus, for example, Central Pattern Gen¬ 
erators (CPGs), i.e. assemblies of small number of 
neurons, capable of providing the necessary rhythm 
of muscular activity even in the absence of ex¬ 
ternal stimuli. These CPGs are believed to play 
an important role in animal locomotion. In these 
CPGs the relevant points to be considered are the 
dynamics of the isolated neurons, e.g. periodic or 
chaotic, the interaction between the oscillators, and 
the way in which information is processed. An 
important aspect is that the resulting spatiotem- 
poral patterns can be analyzed through symmetry 
arguments [Collins & Stewart, 1994], and this 
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allows one to study the different possible behav¬ 
iors, stemming from symmetry-breaking bifurca¬ 
tions, and, thus, the transition between different 
animal gaits has been explained in this way by con¬ 
sidering a model composed out of a ring of coupled 
oscillators [Collins & Stewart, 1994]. Regarding the 
possibility that single neurons are chaotic, some ev¬ 
idences point in this direction [Hayashi & Ishizuka, 
1992], 

In the present work we shall explore further the 
richness of dynamical behaviors that are possible in 
rings of unidirectionally coupled chaotic oscillators, 
considering as case examples the Lorenz and Chua 
systems. We shall perform our study by considering 
a reference state in which the behavior of the cou¬ 
pled systems is chaotic and uniform (synchronized), 
studying then the onset of instability, characterized 
in a Fourier representation by the instability in the 
k — 1 mode. This will yield rotating waves that in 
one case are periodic while in the other chaotic. 


2. Rings of Lorenz Oscillators 

In this case we shall consider rings coupled in such 
a way that the dynamical behavior is defined by, 


ij = a(yj - xj) 

yj = Rxj — yj — Xj Zj > 
Zj =X jyj -bZj 


j = 1,..., N , (1) 


where the coupling enters through Xj, that is de¬ 
fined as xj = Xj- 1 , with x\ = aqy. 

In this situation it was observed [Matfas et al., 
1997a] that the synchronized chaotic state is stable 
if the size of the ring is small enough, e.g. N = 2 
(see Fig. 1), while for a certain critical number, 
iV c = 3 in the case of Lorenz model, an instability 
that destroys the uniform chaotic state occurs, lead¬ 
ing to a rotating chaotic wave (see Fig. 2). We have 
performed studies in which the parameters of the 
system have been varied, with the result that the 
critical size, N c = 3 in all cases. Anyway, so far 
we have explored only the region in which all the 
oscillators are identical. 

A noteworthy aspect of this desynchronization 
transition is that the time scale of the emerging 
rotating wave is, roughly, one order of magnitude 
faster than that of the uncoupled oscillators. This 
instability can be characterized by performing a 
linear stability of the small deviations around the 
synchronized state (see e.g. [Turing, 1952; Heagy 



Fig. 1. Temporal evolution of the variable a: of the two os¬ 
cillators in a ring of N = 2 Lorenz oscillators. The values of 
the parameters are (cr, R, b) = (10, 28, 8/3). 



Fig. 2. Temporal evolution of the variable x of two contigu¬ 
ous oscillators in a ring of N — 3 Lorenz oscillators. Notice 
the time scale of this figure compared to that of Fig. 1. The 
values of the parameters are the same as in Fig. 1. 


et al., 1994]). The time evolution of small differ¬ 
ences around the synchronized state is governed by 
the equation, 

5x = U5x (2) 

where the H matrix is organized in a series of blocks 
corresponding to the uncoupled oscillators plus a 
number of off-diagonal terms arising from coupling. 

However, the structure of this matrix is circu- 
lant, and for this reason one can put these equa¬ 
tions in a more convenient form through the use of 
Discrete Fourier Transform (DFT) [Turing, 1952; 
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Heagy et al, 1994]. As a result, the following equa¬ 
tions are obtained, 


Uq) 


7)( fc ) = C<*> r] (k) (3) 


where in the case of the Lorenz model the structure 
of each block can be cast in the form, 


C (k) 


( —g a 0 \ 

(Re^ — z) —1 —x 

V V x -bj 


( 4 ) 


with efc = exp(i2n k/N) and being k = 0,..., 
(N — 1) the Fourier modes of the system. 

The matrices have time-dependent (chaot¬ 
ically varying) coefficients, and, thus, we have cho¬ 
sen to characterize its stability by determining the 
corresponding Lyapunov spectrum considering the 
infinite-time limit of the real part of the eigenvalues 
of this matrix. This has been done by generaliz¬ 
ing Wolf’s algorithm [Wolf et al., 1985] to the case 
of complex vector spaces, while the different pos¬ 
sible values of k and N have been joined through 
the definition of the reduced wavenumber q = k/N, 
and this yields the function X(q), that represents 
the highest Lyapunov exponent as a function of 
this variable. According to linear stability theory 
[Heagy et al., 1994] the stability of the synchronized 
state will occur whenever the transverse Lyapunov 
exponent is negative. However, if one assumes that 
the dependence of A on q is smooth, the fact that 
the uniform dynamics is chaotic, i.e. that A(0) > 0, 
implies that the uniform chaotic state must be un¬ 
stable for perturbations of some characteristic wave¬ 
length (see also [Bohr et al., 1987] for an analogous 
argument). 

In a more quantitative fashion, it can be shown 
that for the parameters used in this work, and re¬ 
ported in Fig. 3, this crossing occurs for q c ~ 0.37, 
what implies that it occurs already for N = 3. This 
can be confirmed through numerical simulation of 
Eq. (1), as can be seen in Fig. 2. An easily ac¬ 
knowledged point from these results is that when 
the instability occurs the behavior of each oscilla¬ 
tor becomes periodic, and neighboring oscillators 
exhibit a phase difference of 2ix/N. The first aspect, 
i.e. the transition from chaotic to periodic cannot be 
explained in the framework of a linear stability the¬ 
ory, and the observed behavior implies, probably, a 
global bifurcation. 

However, the other aspect can be understood 
by noticing that the instability occurs through a 



q 


Fig. 3. Representation of the highest Lyapunov exponent 
A (q) as a function of q = k/N, i.e. \(q) versus q. The circles 
indicate the highest transverse Lyapunov exponent for a ring 
of N = 2 Lorenz oscillators, whereas the squares indicate the 
same for N = 3. The values of the parameters are the same 
as in Fig. 1. 


symmetric Hopf bifurcation [Collins & Stewart, 
1994]. This bifurcation is allowed because the pres¬ 
ence of the efc terms in Eq. (4) implies that half of 
the Fourier modes are complex conjugate to other 
half. In particular this implies that when a given 
mode crosses the instability threshold there will be 
another mode that also exhibits the same type of 
crossing. Whether these two complex conjugate un¬ 
stable modes are real or complex will depend on 
the structure of the matrix, although in the present 
case, and by resorting to an approximate proce¬ 
dure [Giiemez et al, 1997], we have shown that the 
modes are indeed complex. The result is immediate: 
A Hopf bifurcation occurs, implying the appearance 
of a discrete rotating wave, in which neighboring os¬ 
cillators exhibit the reported phase difference. 

Symmetry is here a very helpful tool as it de¬ 
termines the properties of the different bifurcation 
branches. In the case of unidirectional coupling that 
we are considering here, and that it is the most rel¬ 
evant one in the case of CPGs, symmetry indicates 
that a single branch of rotating waves is obtained 
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[Collins & Stewart, 1994], where a single rotation 
direction is allowed. 


3. Rings of Chua’s Oscillators 

An analogous study to that of the previous section 
has been carried out with rings of Chua’s oscilla¬ 
tors, a well-known paradigm of chaos in electronic 
circuits. In this case the oscillators have been cou¬ 
pled according to the following scheme, where the 
evolution equations for each coupled oscillator are 
reported, 

Xj = a[Vj ~ Xj ~ f(xj )] j 
y =x j -y j + z j > j — 1,..., N . (5) 

Zj=-PVj-'1*j I 


A theoretical study of this situation has been car¬ 
ried out in [Mati'as et al., 1997b], while the pre¬ 
dictions have been confirmed experimentally in 
[Sanchez et al., 1997]. The nonlinear resistor f{x) 
in (5) is given by, 

f{x) = |bx + i(a - 6)[|a; + 1| - |* - 1|]| . (6) 

Driving is introduced through the nonlinear term 
f(x) in (5), such that x k = x k -i for k ^ 1, whereas 
for k = 1 , X\ = xn- 

The same type of linear stability analysis dis¬ 
cussed in the previous section can be applied here. 
In this case the matrix takes the following 

form, 


/ -a[l + f(x)e k ] 


C« = ^ 


1 

0 


a ON 
-1 1 
~P - 7 / 


( 7 ) 


which leads to the representation of the highest 
transverse Lyapunov exponent as a function of the 
reduced wave number, plotted in Fig. 4. In the 
present case it is found that the onset of instability 
occurs at q c = 0 . 21 , which implies that the ring be¬ 
comes unstable when N > 5. This can be seen from 
Fig. 5 that presents results for N = 4, while the 
behavior past the instability is presented in Fig. 6 , 
that shows results for N = 5. The interesting 
feature is now that the behavior of the oscillators 
is not periodic, but chaotic, while neighboring os¬ 
cillators present a phase difference that is approxi¬ 
mately equal to 2rr/N. Thus, the ring can be better 
characterized as exhibiting a rotating chaotic wave. 


Regarding the stability of chaotic synchronization 
as a function of the parameters of the model, 
by increasing a and 7 one sees that the critical 
size N c increases, as shown in Table 1. Thus, 


X(q) 



q 


Fig. 4. Representation of the highest Lyapunov exponent 
A (q) as a function of q = k/N, i.e. \{q) versus q. The cir¬ 
cles indicate the highest transverse Lyapunov exponent for 
a ring of JV = 4 Chua oscillators, whereas the squares indi¬ 
cate the same for N = 5. The values of the parameters are 
(a, P, 7 , a, b) = (10, 14.87, 0.06, -1.27, -0.68). 



Fig. 5. Temporal evolution of the variable x of two contigu¬ 
ous oscillators in a ring of JV = 4 Chua oscillators. The values 
of the parameters are the same as in Fig. 4. 
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Time (t.u.) 

Fig. 6 . Temporal evolution of the variable x of two contigu¬ 
ous oscillators in a ring of IV = 5 Chua oscillators. The values 
of the parameters are the same as in Fig. 4. 


Table 1. Critical number of 
chaotic Chua’s circuits in a 
ring, N c , that supports chaotic 
(uniform) synchronization. 


N c 

(a; 7 ) 

4 

a = 10; 7 < 0.15 

5 

a = 10; 0.15 < 7 < 0.2 

6 

a = 12 ; 7 = 0.2 


for the value 7 = 10 used in all our calculations for 
7 < 0.15 one gets N c = 4, while when 0.15 < 7 < 
0.20 this critical value becomes N c = 5. Consider¬ 
ing higher values of 7 is meaningless, as the system 
becomes periodic. However, one can stabilize the 
IV = 6 ring in the chaotic synchronized state by 
increasing simultaneously a (e.g. for a = 12 and 
7 = 0 . 2 ). 

4. Conclusions 

In the present work we have considered the be¬ 
havior of rings of unidirectionally coupled identical 
chaotic oscillators. In particular, we have consid¬ 
ered the cases of Lorenz and Chua systems. In both 
cases it is found that the uniform chaotic state of 
the ring is unstable to perturbations of some finite 
wavelength, or, equivalently, a finite size N. This 
leads to an instability in the uniform state, that 
is conveniently characterized by performing a Dis¬ 
crete Fourier Transform on the linear stability ma¬ 


trix of the problem. However, the behavior exhib¬ 
ited by these two types of rings is different in that in 
one case (Lorenz system) one obtains discrete peri¬ 
odic rotating waves, while in the second case (Chua 
system) these waves are chaotic. 

These discrete spatiotemporal structures are 
interesting in the context of dynamical systems 
theory, as they reveal the richness of dynamical 
behaviors that one may obtain in coupled arrays 
of chaotic oscillators, beyond synchronization. In 
addition, they can be potentially useful in connec¬ 
tion to CPGs, i.e. rings of coupled neuron models. 
It has been found that these structures can be use¬ 
ful in locomotion, where the different symmetry¬ 
breaking solutions would be responsible for the dif¬ 
ferent gaits that the animal exhibits [Collins & 
Stewart, 1994]. In particular, in this context it is 
useful to notice that in the case of rings of Lorenz 
systems the pattern that emerges is periodic, al¬ 
though the dynamics of the uncoupled oscillators 
was chaotic. It is also interesting to notice the time 
scale of the emerging rotating wave: It is, at least, 
one order of magnitude faster than the uncoupled 
oscillators, and this could be relevant in a neuronal 
context. As well, one should bear in mind that the 
brain is able to perform various tasks in a short 
time, although the neurons in which these tasks 
base are relatively slow. 
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I report on the experimental observation of excitation and suppression of chaos through time 
dependent perturbations in the dynamical variable of a glow discharge. The interaction of the 
external signal with the dynamical system is explained in terms of the ID map associated to the 
glow discharge. Numerical simulations are also performed with the logistic map. The proposed 
mechanism of exciting and/or suppressing chaos is in accordance with the OGY method of 
controlling chaos. 


1. Introduction 

Very often a dynamical system may be modeled by a 
discrete time evolution equation that can be cast in 
the form Xj+i = /(x,, fi). This is a generic example 
of a unidimensional (ID) map, where x, represents 
the dynamical variable and /i is a control parame¬ 
ter. It is well known that by varying p, the system 
may present bifurcations (i.e. qualitative changes in 
its dynamical behavior) leading to a very complex 
evolution that we assume to be chaotic. The prop¬ 
erties of this road to chaos have been extensively 
studied in the literature for several ID mappings, 
as for instance the logistic map: Xj+i = px,( 1 — Xj). 

In this work, besides the logistic map, an 
experimental system is analyzed: the glow dis¬ 
charge. Experimental systems are dissipative, thus 
their phase-space volumes contract in all directions. 
One direction has the slowest convergence and de¬ 
fines a line along which the universality theory ap¬ 
plies [Collet & Eckmann, 1980] and therefore their 
chaotic dynamics can be modeled by ID maps. 

In the example of the glow discharge the dy¬ 
namical variable is the electric current flowing 
through the discharge and the control parameter is 
the DC voltage feeding the discharge. Under suit¬ 
able conditions the current of the discharge presents 
a self-generated periodic oscillation that, by chang¬ 


ing the applied voltage, may turn into a chaotic 
oscillation. The associated bifurcation sequences 
eventually reaching chaos are characterized else¬ 
where [Braun et al., 1987, 1992, 1994, 1995]. In 
this paper suppression and excitation of chaos in 
the glow discharge are analyzed by using time- 
dependent perturbations on the current of the dis¬ 
charge. The mechanism of control is explained in 
terms of the ID map associated to the discharge 
and complemented with the example of the logistic 
map. 

This paper is organized as follows: in Sec. 2, 
I describe the experimental apparatus. The mea¬ 
surements and their analysis are presented in Sec. 3. 
Conclusions are presented in Sec. 4. 

2. Experimental 

The experimental arrangement is shown schemati¬ 
cally in Fig. 1. The main setup consists of a dis¬ 
charge cell connected to an adjustable DC power 
supply (Ortec, model 456) and a ballast resistor of 
1 Mfl in a series circuit. The discharge cell is as¬ 
sembled with a glass tube 24 mm in diameter and 
has brass electrodes which are ~ 1 cm apart. The 
cathode has a diameter of 12 mm while the anode 
has a diameter of 19.2 mm. Inside the cell there 
is an argon pressure of ~ 2 mbar. The cathode is 
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Fig. 1. Experimental setup. Vs, power supply; Rs = 1 MSI; 
Ri = 2 kfi; DC, discharge cell; DB, digitizing board; FG, 
function generator. 


connected to the voltage source and the anode is 
grounded through a 2 M2 resistor which is used to 
measure the current. The current was monitored 
by a digitizing board (Sonix STR825) in a personal 
microcomputer. The magnitude of this current is of 
the order of some fiA when the power supply fur¬ 
nishes a voltage that ranges between 300 and 500 V. 

Depending on the value of this voltage, the cur¬ 
rent through the discharge shows periodic or chaotic 
oscillations. In order to excite or suppress chaos in 
the glow discharge the current is perturbed by an 
external signal. The perturbation arrangement is 
shown in Fig. 1 by the drawing in gray. The per¬ 
turbation signal was applied to the discharge by a 
capacitive coupling. The output from a variable 
frequency function generator (Tektronix FG 504) 
was connected to a copper tube of 20 mm in length 
which just fitted the outside of the glass discharge 
cell. In this way the electric field distribution inside 
the discharge cell is perturbed by the signal gener¬ 
ator. I adopted the following criterion to define a 
small perturbation. With the discharge out of oper¬ 
ation I measured the voltage induced by the exter¬ 
nal signal between electrodes. This voltage must be 
less than 5% of the voltage furnished by the power 
supply when the discharge is in operation. As will 
be explained later, the perturbation signal may be 
periodic or not. 

3. Results and Analysis 

An example of the chaotic evolution of the current 
in the discharge is displayed at the left of Fig. 2, 
whereas at the right the corresponding next am¬ 
plitude map is shown. This map is obtained by 
sampling only the amplitudes (A;) of the current 
oscillations and arranging them in a plot A* x Ai + \. 
The result is a noticeable ID map. 



(a) 



(b) 


Fig. 2. (a) Chaotic current in the glow discharge. Two suc¬ 

cessive amplitudes (A, and A,+i) are identified, (b) Next 
amplitude map corresponding to (a). 


By perturbing the electrical discharge with an 
external periodic signal, the dynamics changes dras¬ 
tically. This can be seen in Figs. 3 and 4. Both 
diagrams are obtained digitizing only the ampli¬ 
tudes of the oscillations. The continuous lines 
in this diagram characterize a periodic evolution 
whereas a blurred collection of points indicates a 
chaotic evolution. Also, both diagrams show a dras¬ 
tic change in dynamics when the perturbation is 
switched on and off. While the perturbation is ac¬ 
tive, its amplitude is slowly increased linearly. 

In Fig. 3 the initially periodic oscillations (pe¬ 
riod one) on the current make an abrupt transi¬ 
tion to chaos with a very small external perturba¬ 
tion demonstrating the excitation of chaos. As the 
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Fig. 3. Excitation of chaos in the glow discharge. The two 
arrows indicate where the perturbation is turned, respec¬ 
tively, on and off. 
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Fig. 4. Suppression of chaos in the glow discharge. 


amplitude of the perturbation is increased, a se¬ 
quence of oscillations with period .. .-8-7-6-5-... is 
clearly observed between chaotic regions. When the 
external signal is switched off, the current returns 
to its initial state. 

The suppression of chaos is shown in Fig. 4. 
The initially chaotic behavior of the current changes 
to an oscillation of period four as the amplitude 
of the external signal is increased. The oscillation 
undergoes a reverse period doubling, resulting in 
a period-two oscillation at a higher amplitude of 
the perturbation. The current returns to its ini¬ 
tial noisy appearance when the external oscillator 
is switched off. 

Similar results are obtained for different exter¬ 
nal periodic signals: a sine wave, a triangular wave, 
and a square wave. All of them have a frequency 
very close to the one of the auto-oscillations of the 
current. Also, a nonperiodic signal is effective in 


exciting and/or suppressing chaos. For example, 
instead of using the function generator in the free 
run mode, it may be triggered such that the out¬ 
put of the generator is always synchronized with 
the amplitudes of the oscillation. For this purpose, 
the trigger input was generated by injecting the dis¬ 
charge current oscillations in a peak detector (differ¬ 
entiator plus comparator in Fig. 1). In this way, the 
perturbation has always the same periodicity of the 
oscillations in the discharge current. For instance, if 
the current is chaotic then the perturbation signal is 
also chaotic, but as the the perturbation increases 
and achieves suppression of chaos then the signal 
becomes periodic. 

I used an external perturbation in the form of 
“delta” pulses synchronized with the amplitudes of 
the current oscillation to excite and/or suppress 
chaos. In this way the perturbation systemati¬ 
cally changes the variable of the next amplitude 
map. The latter observation suggests that the in¬ 
teraction of the external signal with the dynami¬ 
cal system (discharge) can best be understood in 
terms of the ID map associated to the system as 
the next amplitude map. I empirically propose 
that this interaction takes then the form: Xj+i = 
f(xi, n) Inconstant — £i represents the external 
perturbation; its amplitude increases linearly with 
the successive iterations. This linear growth better 
fits the experimental results presented previously. 

To check this proposition, I have made some nu¬ 
merical simulations using the logistic map: Xi+i = 
fxxi( 1 — Xi) + < j £ — ei (see Fig. 5). Concerning 
the glow discharge, this map may be considered as 
“equivalent” to the ID map ruling the dynamics 
such as the next amplitude map or a Poincare re¬ 
turn map that can be obtained from the phase 



n (iterations) 

Fig. 5. Simulation for the logistic map with p = 3.75. 
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space flow through a Poincare section. I considered 
Si = si, s being constant and i corresponding to the 
iteration. A Gaussian noise <r£ (a = amplitude = 
0.002, (£) = 0; (£ 2 ) = 1) was introduced in another 
term to represent external noise, always present in 
experimental conditions. The noise does not dis¬ 
turb the proposed control mechanism. It may be 
concluded that the dynamical state of the system 
with n = constant varies according to the ampli¬ 
tude of the perturbation £j. In fact, £{ works as a 
“new” control parameter and, changing linearly its 
amplitude, it results a similar bifurcation sequence 
as that observed in the nonperturbed system. 

Related procedures of perturbing the vari¬ 
ables in ID discrete time systems (maps) have 
been studied by Giiemez and Matfas [1993] and 
Parthasarathy and Sinha [1995]. The present work, 
besides demonstrating new features of such ideas in 
ID maps, extends them to continuous time systems 
(flows) like the glow discharge. The mechanism of 
exciting and suppressing chaos outlined in this pa¬ 
per has also been applied to a numerical Rossler 
system. It turned out that this control procedure is 
also effective in a numerical continuous time system 
[Braun, 1997]. 

As a final remark, I want to stress that the 
proposed mechanism of exciting and/or suppress¬ 
ing chaos is in accordance with the OGY method 
of controlling chaos [Ott et al., 1990; Shinbrot 
et al., 1993], which consists of stabilizing an unsta¬ 
ble periodic orbit present in the system under con¬ 
trol. The mechanism proposed in this work does not 
change the control parameter, it only displaces up- 
or downwards the entire function f(xi, //) of the ID 
map by an amount Si. The direction of the displace¬ 
ment depends on the sign of e. Possible period-n 
orbits of the map (n = 1, 2, 3,...) are given by the 
intersection of with the diagonal y = x. They 
are stable or unstable whether f' 71 ' 1 (x,) is respec¬ 
tively greater or less than one. By displacing f( n \ 
the stability of these periodic points is changed and 
thus chaos can be controlled by stabilizing an unsta¬ 
ble periodic orbit, what is the essential idea of the 
OGY method. In the same way a stable periodic 
orbit can be turned into an unstable one. 


4. Conclusions 

I demonstrated that it is possible to suppress and/or 
excite chaos through time dependent perturbations 
in the variable of a dynamical system. I investigated 
experimentally the glow discharge and numerically 
the logistic map. The proposed control mechanism 
is easily understood in terms of the ID map asso¬ 
ciated to the dynamical system. Essentially it con¬ 
sists in changing the stability of the periodic points 
of the map. 
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We introduce the idea of the fuzzy control of chaos: we show how fuzzy logic can be applied to 
the control of chaos, and provide an example of fuzzy control used to control chaos in Chua’s 
circuit. 


1. Introduction 

Chaos control exploits the sensitivity to initial con¬ 
ditions and to perturbations that is inherent in 
chaos as a means to stabilize unstable periodic or¬ 
bits within a chaotic attractor. The control can 
operate by altering system variables or system 
parameters, and either by discrete corrections or by 
continuous feedback. Many methods of chaos con¬ 
trol have been derived and tested [Chen & Dong, 
1993; Lindner & Ditto, 1995; Ogorzalek, 1993]. 
Why then consider fuzzy control of chaos? 

A fuzzy controller works by controlling a con¬ 
ventional control method. We propose that fuzzy 
control can become useful together with one of these 
other methods — as an extra layer of control — in 
order to improve the effectiveness of the control in 
terms of the size of the region over which control is 
possible, the robustness to noise, and the ability to 
control long period orbits. 

In this paper, we put forward the idea of fuzzy 
control of chaos, and we provide an example show¬ 
ing how a fuzzy controller applying occasional pro¬ 
portional feedback to one of the system parameters 
can control chaos in Chua’s circuit. 


2. Fuzzy Control 

Fuzzy control [Driankov et al, 1993; Terano el al., 
1994] is based on the theory of fuzzy sets and fuzzy 
logic [Yager & Zadeh, 1991; Bezdek, 1993]. The 
principle behind the technique is that imprecise 
data can be classified into sets having fuzzy rather 
than sharp boundaries, which can be manipulated 
to provide a framework for approximate reason¬ 
ing in the face of imprecise and uncertain infor¬ 
mation. Given a datum, x, a fuzzy set A is said 
to contain x with a degree of membership 
where ha{x) can take any value in the domain 
[0,1]. Fuzzy sets are often given descriptive names 
(called linguistic variables) such as FAST; the mem¬ 
bership function //fast/*) is then used to reflect 
the similarity between values of x and a contex¬ 
tual meaning of FAST. For example, if x repre¬ 
sents the speed of a car in kilometres per hour, 
and FAST is to be used to classify cars travelling 
fast, then FAST might have a membership func¬ 
tion equal to zero for speeds below 90 km/h and 
equal to one for speeds above 130 km/h, with a 
curve joining these two extremes for speeds between 
these values. The degree of truth of the statement 
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the car is travelling fast is then evaluated by read¬ 
ing off the value of the membership function corre¬ 
sponding to the car’s speed. 

Logical operations on fuzzy sets require an ex¬ 
tension of the rules of classical logic. The three 
fundamental Boolean logic operations, intersection, 
union, and complement, have fuzzy counterparts 
defined by extension of the rules of Boolean logic. A 
fuzzy expert system uses a set of membership func¬ 
tions and fuzzy logic rules to reason about data. 
The rules are of the form “if x is FAST and y is 
SLOW then 2 is MEDIUM”, where x and y are in¬ 
put variables, 2 is an output variable, and SLOW, 
MEDIUM, and FAST are linguistic variables. The 
set of rules in a fuzzy expert system is known as the 
rule base, and together with the data base of input 
and output membership functions it comprises the 
knowledge base of the system. 

A fuzzy expert system functions in four steps. 
The first is fuzzification, during which the member¬ 
ship functions defined on the input variables are ap¬ 
plied to their actual values, to determine the degree 
of truth for each rule premise. Next under inference, 
the truth value for the premise of each rule is com¬ 
puted, and applied to the conclusion part of each 
rule. This results in one fuzzy set to be assigned to 
each output variable for each rule. In composition, 
all of the fuzzy sets assigned to each output vari¬ 
able are combined together to form a single fuzzy 
set for each output variable. Finally comes defuzzi¬ 
fication, which converts the fuzzy output set to a 
crisp (nonfuzzy) number. 

A fuzzy controller may then be designed using 
a fuzzy expert system to perform fuzzy logic opera¬ 
tions on fuzzy sets representing linguistic variables 
in a qualitative set of control rules (see Fig. 1). 

As a simple metaphor of fuzzy control in prac¬ 
tice, consider the experience of balancing a stick 



Fig. 1. Fuzzy logic controller block diagram. 


vertically on the palm of one’s hand. The equations 
of motion for the stick (a pendulum at its unstable 
fixed point) are well-known, but we do not inte¬ 
grate these equations in order to balance the stick. 
Rather, we stare at the top of the stick and carry 
out a type of fuzzy control to keep the stick in the 
air: We move our hand slowly when the stick leans 
by a small angle, and fast when it leans by a larger 
angle. Our ability to balance the stick despite the 
imprecision of our knowledge of the system is at the 
heart of fuzzy control. 


3. Techniques for Fuzzy Chaos 
Control 

To control a system necessitates perturbing it. 
Whether to perturb the system via variables or pa¬ 
rameters depends on which are more readily acces¬ 
sible to be changed, which in turn depends on what 
type of system is to be controlled — electronic, me¬ 
chanical, optical, chemical, biological, etc. Whether 
to perturb continuously or discretely is a question 
of intrusiveness — it is less intrusive to the system, 
and less expensive to the controller, to perturb dis¬ 
cretely. Only when discrete control is not effective 
might continuous control be considered. 

Ott, Grebogi and Yorke [Ott et al, 1990] 
invented a method of applying small feedback 
perturbations to an accessible system parameter in 
order to control chaos. The OGY method uses the 
dynamics of the linearized map around the orbit 
one wishes to control. Using the OGY method, 
one can pick any unstable periodic orbit that exists 
within the attractor and stabilize it. The control is 
imposed when the orbit crosses a Poincare section 
constructed close to the desired unstable periodic 
orbit. Since the perturbation applied is small, it is 
supposed that the unstable periodic orbit is unaf¬ 
fected by the control. 

Occasional proportional feedback [Hunt, 1991; 
Lindner & Ditto, 1995] is a variant of the original 
OGY chaos control method. Instead of using the 
unstable manifold of the attractor to compute cor¬ 
rections, it uses one of the dynamical variables, in 
a type of one-dimensional OGY method. This feed¬ 
back could be applied continuously or discretely in 
time; in occasional proportional feedback it is ap¬ 
plied discretely. Occasional proportional feedback 
exploits the strongly dissipative nature of the flows 
often encountered, enabling one to control them 
with a one-dimensional map. The method is easy 
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to implement, and in many cases one can stabi¬ 
lize high period unstable orbits by using multiple 
corrections per period. It is a suitable method on 
which to base a fuzzy logic technique for the control 
of chaos, since it requires no knowledge of a system 
model, but merely an accessible system parameter. 


4. An Example: Fuzzy Control 
of Chaos in Chua’s Circuit 

Chua’s circuit [Matsumoto, 1984; Kennedy, 1993] 
exhibits chaotic behavior that has been extensively 
studied, and whose dynamics is well known [Madan, 
1993]. Recently, occasional proportional feedback 
has been used to control the circuit [Johnson et al. 
1993]. The control used an electronic circuit to sam¬ 
ple the peaks of the voltage across the negative re¬ 
sistance and if it fell within a window, centred about 
a by a set-point value, modified the slope of the 
negative resistance by an amount proportional to 
the difference between the set point and the peak 


NEGATIVE NEGATIVE POSITIVE POSITIVE 

BIG SMALL ZERO SMALL BIG 



Fig. 2. Membership functions of the input and output 
variables e, Ae, and A a. 

value. The nonlinear nature of this system and the 
heuristic approach used to find the best set of pa¬ 
rameters to take the system to a given periodic or¬ 
bit suggest that a fuzzy controller that can include 
knowledge rules to achieve periodic orbits may pro¬ 
vide significant gains over occasional proportional 
feedback alone. 

We have implemented a fuzzy controller to con¬ 
trol the nonlinearity of the nonlinear element (a 
three segment nonlinear resistance) within Chua’s 
circuit. The block diagram of the controller is 



Fig. 3. The whole controller and control system in the form of a block diagram, including the fuzzy controller, the peak 
detector, the window comparator, and the Chua’s circuit system being controlled. 


Table 1. Quantification levels and membership functions. 


Error, e 

-1 

-0.75 

m 

-0.25 

0 

0.25 
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Change in error, Ae 

-1 
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Linguistic Variables Membership Functions 


Positive Big, PB 

0 

0 

0 

0 

0 

0 

0 
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Positive Small, PS 

0 

0 

0 

0 

0 
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1 

0.5 
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0 

0 

0 

0.5 

1 

0.5 

0 

0 

0 

Negative Small, NS 

0 

0.5 
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0.5 

0 

0 

0 

0 

0 

Negative Big, NB 

1 

0.5 

0 

0 

0 

0 

0 

0 

0 
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Table 2. Rule table for the linguistic vari¬ 
ables in Table 1. 


\ e 
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AZ 
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AZ 

NS 

NS 

AZ 

AZ 

PS 


AZ 

NS 

AZ 
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PS 


PS 

AZ 

AZ 

PS 

PB 


PB 

AZ 

PS 

PS 

PB 



shown in Fig. 1. It consists of four blocks: knowl¬ 
edge base, fuzzification, inference and defuzzifica¬ 
tion. The knowledge base is composed of a data 
base and a rule base. The data base consists of the 
input and output membership functions (Fig. 2). 
It provides the basis for the fuzzification, defuzzi¬ 
fication and inference mechanisms. The rule base 
is made up of a set of linguistic rules mapping in¬ 
puts to control actions. Fuzzification converts the 
input signals e and Ae into fuzzified signals with 
membership values assigned to linguistic sets. The 
inference mechanisms operate on each rule, apply¬ 
ing fuzzy operations on the antecedents and by 
compositional inference methods derives the con¬ 
sequents. Finally, defuzzification converts the fuzzy 
outputs to control signals, which in our case control 
the slope of the negative resistance Aa in Chua’s 



) 1 ...■- 1 — 

0.0 100.0 
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Fig. 4. The fuzzy controller stabilizes a previously unsta¬ 
ble period-1 orbit. The control is switched on at time 20. 
The lower trace shows the correction pulses applied by the 
controller. 

circuit (Fig. 3). The fuzzification maps the error e, 
and the change in the error Ae, to labels of fuzzy 
sets. Scaling and quantification operations are ap¬ 
plied to the inputs. Table 1 shows the quantified 
levels and the linguistic labels used for inputs and 
output. The knowledge rules (Table 2) are repre¬ 
sented as control statements such as “if e is NEG¬ 
ATIVE BIG and Ae is NEGATIVE SMALL then 
Aa is NEGATIVE BIG”. 



0.0 100.0 


(a) (b) 

Fig. 5. Trajectory traces show higher period orbits stabilized by the controller. As before, the lower trace shows the correction 
pulses applied by the fuzzy control. 
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The normalized equations representing the cir¬ 
cuit are 

x = a{y - x - f(x )), 
y = x — y + z, (1) 

i = -Py, 

where f(x) =bx + |(a — 6)(|m + 1| — \x — 1|) repre¬ 
sents the nonlinear element of the circuit. Changes 
in the negative resistance were made by changing a 
by an amount 

A a = Fuzzy Controller Output x Gain x a. (2) 

We have performed numerical simulations, both 
in C and in Simulink, of Chua’s circuit controlled 
by the fuzzy logic controller. Figure 3 shows the 
whole control system in the form of a block dia¬ 
gram, including Chua’s circuit, the fuzzy controller, 
the peak detector, and the window comparator. 
Figure 4 gives a sample output of the fuzzy con¬ 
troller stabilizing an unstable period-1 orbit by ap¬ 
plying a single correction pulse per cycle of oscilla¬ 
tion. By changing the control parameters we can 
stabilize orbits of different periods. In Fig. 5 we 
illustrate more complex higher period orbits stabi¬ 
lized by the controller. One can tune the fuzzy con¬ 
trol over the circuit to achieve the type of response 
required in a given situation by modifying some or 
all of the rules in the knowledge base of the system. 

Of course, in the case of Chua’s circuit the sys¬ 
tem equations are available and fuzzy logic is thus 
not necessary for control, but this simple example 
permits us to see the possibilities that fuzzy control 
provides, by allowing a nonlinear gain implemented 
in the form of knowledge based rules. 

5. Conclusions 

We have introduced the idea of using fuzzy logic 
for the control of chaos. Fuzzy logic controllers are 
commonly used to control systems whose dynamics 


is complex and unknown, but for expositional clar¬ 
ity here we have given an example of its use with a 
well-studied chaotic system. We have shown that it 
is possible to control chaos in Chua’s circuit using 
fuzzy control. Further work is necessary to quantify 
the effectiveness of fuzzy control of chaos compared 
with alternative methods, to identify ways in which 
to systematically build the knowledge base for fuzzy 
control of a particular chaotic system, and to apply 
the fuzzy controller to future chaotic systems. 
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